mk_grid_c2f

MK_GRID_C2F - calculate a coarse2fine mapping for grid coarse models.
 C2F = MK_GRID_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 (vox-based) model.
 Uses CONVHULLN to calculate the volume defined by a set of intersection
 points between individual tet and vox elements.

 C2F = MK_GRID_C2F(FMDL,RMDL,OPT) allows specifying options.
 
 Inputs:
   FMDL - an EIDORS (tet-based) forward model
   RMDL - a grid model, as returned by MK_GRID_MODEL
   OPT  - an option structure with the following fields and defaults:
      .do_not_scale  - set to true to prevent scaling the models to unit
                       cube before any calculations, including thresholds.
                       Default: false
      .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 a point is lying inside a triangle.
                       Default: eps
      .save_memory   - modifies function behavior to decrease memory 
                       footprint by increasing the number of iterations;
                       useful for large problems. Must be an integer 
                       between 0 and 3
                          0  -  calculate all at once (default)
                          1  -  calculate one xy voxel plane at a time
                          2  -  calculate one y voxel row at a time
                          3  -  calculate each voxel separately
                       NOTE: read below about persistent usage

 MK_GRID_C2F('save_memory',N) sets the 'save memory' option for all future
 calls (persistent variable). N will be used when opt.save_memory is not
 specified.

 [C2F M] = MK_GRID_C2F(...) also returns a struct with useful fields
 characterising the vox model

 Set eidors_msg 'log level' < 2 to supress output to command line.

 Examples:
     fmdl = ng_mk_cyl_models([2,2,.2],[],[]);
     rmdl = mk_grid_model([],-2:2,-2:2,0:2);
     c2f  = mk_grid_c2f(fmdl,rmdl);
     h = show_fem(fmdl); set(h,'LineWidth',0.1)
     hold on
     h = show_fem(rmdl); set(h,'EdgeColor','b','LineWidth',2);
     hold off

 See also MK_GRID_MODEL, FIND_EDGE2EDGE_INTERSECTIONS, CONVHULLN
     MK_TET_C2F, MK_APPROX_C2F, POINT_IN_TRIANGLE, EIDORS_MSG

0001 function [c2f, m] = mk_grid_c2f(fmdl, rmdl, opt)
0002 %MK_GRID_C2F - calculate a coarse2fine mapping for grid coarse models.
0003 % C2F = MK_GRID_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 (vox-based) model.
0006 % Uses CONVHULLN to calculate the volume defined by a set of intersection
0007 % points between individual tet and vox elements.
0008 %
0009 % C2F = MK_GRID_C2F(FMDL,RMDL,OPT) allows specifying options.
0010 %
0011 % Inputs:
0012 %   FMDL - an EIDORS (tet-based) forward model
0013 %   RMDL - a grid model, as returned by MK_GRID_MODEL
0014 %   OPT  - an option structure with the following fields and defaults:
0015 %      .do_not_scale  - set to true to prevent scaling the models to unit
0016 %                       cube before any calculations, including thresholds.
0017 %                       Default: false
0018 %      .tol_node2tet  - tolerance for determinant <= 0 in testing for
0019 %                       points inside tets. Default: eps
0020 %      .tol_edge2edge - maximum distance between "intersecting" edges
0021 %                       Default: 6*sqrt(3)*eps(a), where a is
0022 %                       min(max(abs(fmdl.nodes(:))),max(abs(rmdl.nodes(:)))
0023 %      .tol_edge2tri  - minimum value of a barycentric coordinate to
0024 %                       decide a point is lying inside a triangle.
0025 %                       Default: eps
0026 %      .save_memory   - modifies function behavior to decrease memory
0027 %                       footprint by increasing the number of iterations;
0028 %                       useful for large problems. Must be an integer
0029 %                       between 0 and 3
0030 %                          0  -  calculate all at once (default)
0031 %                          1  -  calculate one xy voxel plane at a time
0032 %                          2  -  calculate one y voxel row at a time
0033 %                          3  -  calculate each voxel separately
0034 %                       NOTE: read below about persistent usage
0035 %
0036 % MK_GRID_C2F('save_memory',N) sets the 'save memory' option for all future
0037 % calls (persistent variable). N will be used when opt.save_memory is not
0038 % specified.
0039 %
0040 % [C2F M] = MK_GRID_C2F(...) also returns a struct with useful fields
0041 % characterising the vox model
0042 %
0043 % Set eidors_msg 'log level' < 2 to supress output to command line.
0044 %
0045 % Examples:
0046 %     fmdl = ng_mk_cyl_models([2,2,.2],[],[]);
0047 %     rmdl = mk_grid_model([],-2:2,-2:2,0:2);
0048 %     c2f  = mk_grid_c2f(fmdl,rmdl);
0049 %     h = show_fem(fmdl); set(h,'LineWidth',0.1)
0050 %     hold on
0051 %     h = show_fem(rmdl); set(h,'EdgeColor','b','LineWidth',2);
0052 %     hold off
0053 %
0054 % See also MK_GRID_MODEL, FIND_EDGE2EDGE_INTERSECTIONS, CONVHULLN
0055 %     MK_TET_C2F, MK_APPROX_C2F, POINT_IN_TRIANGLE, EIDORS_MSG
0056 
0057 % (C) 2015 Bartlomiej Grychtol - all rights reserved by Swisstom AG
0058 % License: GPL version 2 or 3
0059 % $Id: mk_grid_c2f.m 7040 2024-11-30 00:02:08Z bgrychtol $
0060 
0061 persistent save_memory;
0062 
0063 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'), do_unit_test; return; end
0064 if ischar(fmdl) && strcmp(fmdl,'PROFILE'), do_small_test; return; end
0065 if ischar(fmdl) && strcmp(fmdl,'save_memory')
0066    if nargin ~= 2, error('Expected a value for save_memory option'); end
0067    if ischar(rmdl), rmdl = str2double(rmdl); end
0068    save_memory = rmdl;
0069    return
0070 end
0071 if nargin < 3
0072    opt = struct;
0073 end
0074 if ~isempty(save_memory)
0075    try 
0076       opt.save_memory; % command line options have precedence
0077    catch
0078       opt.save_memory = save_memory;
0079    end
0080 end
0081 
0082 [fmdl,rmdl] = eidors_cache(@center_scale_models,{fmdl,rmdl, opt});
0083 
0084 opt = eidors_cache(@parse_opts,{fmdl,rmdl, opt});
0085 
0086 copt.cache_obj = {fmdl.nodes,
0087                   fmdl.elems,
0088                   rmdl.nodes,
0089                   rmdl.elems,
0090                   rmfield(opt,'save_memory')};
0091                
0092 copt.fstr = 'mk_grid_c2f';
0093 
0094 [c2f, m] = eidors_cache(@do_mk_grid_c2f,{fmdl,rmdl,opt},copt);
0095 
0096 
0097 
0098 %-------------------------------------------------------------------------%
0099 % Wrapper providing the save_memory functionality
0100 function [c2f, m]= do_mk_grid_c2f(fmdl0,rmdl0,opt0)
0101     eidors_msg('@@@',2);
0102     
0103     m = [];
0104     c2f = sparse(opt0.nTet,opt0.nVox);
0105    
0106     progress_msg('Prepare tet model...');
0107     fmdl0 = prepare_fmdl(fmdl0);
0108     progress_msg(Inf);
0109     
0110     if opt0.save_memory == 0
0111        relem_idx = ones(opt0.nVox,1);
0112        [fmdl, rmdl, opt, felem_idx] = crop_models(fmdl0,rmdl0,opt0,relem_idx);
0113        if any(felem_idx)
0114           [tmp, m] = separable_calculations(fmdl,rmdl0,opt);
0115           c2f = combine_c2f(c2f, tmp,felem_idx,relem_idx);
0116        end
0117     elseif opt0.save_memory == 10
0118 % The idea here is to calculate separately on each fraction
0119 % of the model. It should then be possible to run in a parfor
0120 % loop, ... but matlab has subtle bugs
0121        logmsg(' Saving memory mode level %d\n',opt0.save_memory);
0122        n_elems = num_elems(fmdl0); step = 5e4;
0123        max_iter = floor(n_elems/step);
0124        pmopt.final_msg = '';
0125        ctr = interp_mesh(fmdl0);
0126        [~,xidx] = sort(ctr(:,1)); % sort by x - real solution is travelling salesman
0127        [xl] = ndgrid(opt0.xvec(1:end-1),opt0.yvec(1:end-1),opt0.zvec(1:end-1)); xl=xl(:);
0128        [xu] = ndgrid(opt0.xvec(2:end-0),opt0.yvec(2:end-0),opt0.zvec(2:end-0)); xu=xu(:);
0129        progress_msg('Progress:',0,max_iter+1,pmopt);
0130        progress = 0;
0131        for k = 0:max_iter;
0132           eidx = true(n_elems,1);
0133           eidx( xidx((k*step+1):min((k+1)*step,n_elems)) ) = false;
0134           fmdl = fmdl0; fmdl.elems(eidx,:) = [];
0135                         fmdl.edge2elem(:,eidx)= [];
0136                         fmdl.node2elem(:,eidx)= [];
0137                         fmdl.elem2face(eidx,:)= [];
0138           opt = opt0;   opt.nTet = sum(eidx==0);
0139 if 1
0140           xvals = fmdl.nodes(fmdl.elems(:),1);
0141           ridx = true(size(xu));
0142           ridx(xl>max(xvals)) = false;
0143           ridx(xu<min(xvals)) = false;
0144           rmdl = rmdl0; 
0145              idx = rmdl0.coarse2fine * ridx > 0;
0146              rmdl.elems = rmdl0.elems(idx,:);
0147              [node_idx, ~,Nn] = unique(rmdl.elems(:));
0148              rmdl.elems = reshape(Nn,size(rmdl.elems));
0149              rmdl.nodes = rmdl0.nodes(node_idx,:);
0150           opt.xvec = unique(rmdl.nodes( rmdl.elems(:), 1));
0151           opt.Xsz = length(opt.xvec)-1;
0152           opt.ystep = opt.xstep*opt.Xsz+1;
0153           opt.zstep = opt.ystep*(opt.Ysz+1);
0154           opt.nVox = opt.Xsz*opt.Ysz*opt.Zsz;
0155 
0156    if 0 % Display the cut model for debugging
0157           fmdl.boundary = find_boundary(fmdl);
0158           rmdl.boundary = find_boundary(rmdl);
0159           hh=show_fem(rmdl); set(hh,'EdgeColor',[0,0,1]);
0160               hold on; show_fem(fmdl); hold off; keyboard
0161    end
0162 else    % Test code: don't cut the rmdl - slower but correct
0163           rmdl = rmdl0; ridx= true(opt.nVox,1);
0164 end
0165           if isempty(rmdl.elems); continue; end % don't bother if rmdl outside
0166 
0167           eidors_msg('log_level',eidors_msg('log_level')-2);
0168           c2f(~eidx,ridx) = separable_calculations(fmdl,rmdl,opt);
0169           eidors_msg('log_level',eidors_msg('log_level')+2);
0170           progress_msg(k+1, max_iter+1);
0171        end
0172        progress_msg(Inf);
0173     else % save_memory > 0
0174        logmsg(' Saving memory mode level %d\n',opt0.save_memory);
0175        max_iter = opt0.Zsz;
0176        if opt0.save_memory >= 2, max_iter = max_iter * opt0.Ysz; end
0177        if opt0.save_memory == 3, max_iter = max_iter * opt0.Xsz; end
0178        pmopt.final_msg = '';
0179        progress_msg('Progress:',0,max_iter,pmopt);
0180        progress = 0;
0181        opt = opt0;
0182        m = prepare_vox_mdl(rmdl0,opt0);
0183        for z = 1:opt0.Zsz
0184           opt.Zsz = 1;
0185           opt.zvec = opt0.zvec(z:z+1);
0186           opt.xplane = opt0.xplane(z,:);
0187           opt.yplane = opt0.yplane(z,:);
0188           relem_idx = false(opt0.Xsz*opt0.Ysz*opt0.Zsz,1);
0189           relem_idx((z-1)*opt0.Xsz*opt0.Ysz + (1:opt0.Xsz*opt0.Ysz)) = true;
0190           if opt0.save_memory > 1
0191              for y = 1:opt0.Ysz
0192                 opt.Ysz = 1;
0193                 opt.yvec = opt0.yvec(y:y+1);
0194                 opt.zplane = opt0.zplane(y,:);
0195                 opt.xplane = opt0.xplane(z,y);
0196                 opt.zstep = opt.ystep*(opt.Ysz+1);
0197                 relem_idx_y = false(opt0.Xsz*opt0.Ysz,1);
0198                 relem_idx_y((y-1)*opt0.Xsz + (1:opt0.Xsz)) = true;
0199                 relem_idx_y = relem_idx & repmat(relem_idx_y,opt0.Zsz,1);
0200                 if opt0.save_memory > 2
0201                    for x = 1:opt0.Xsz
0202                       opt.Xsz = 1;
0203                       opt.xvec = opt0.xvec(x:x+1);
0204                       opt.yplane = opt0.yplane(z,x);
0205                       opt.zplane = opt0.zplane(y,x);
0206                       opt.ystep = opt.xstep*opt.Xsz+1;
0207                       opt.zstep = opt.ystep*(opt.Ysz+1);
0208                       relem_idx_x = false(opt0.Xsz,1);
0209                       relem_idx_x(x) = true;
0210                       relem_idx_x = relem_idx_y & repmat(relem_idx_x,opt0.Ysz*opt0.Zsz,1);
0211                       [fmdl, rmdl, opt, felem_idx] = crop_models(fmdl0,rmdl0,opt,relem_idx_x);
0212                       if any(felem_idx)
0213                          eidors_msg('log_level',eidors_msg('log_level')-2);
0214                          tmp = separable_calculations(fmdl,rmdl,opt);
0215                          eidors_msg('log_level',eidors_msg('log_level')+2);
0216                          c2f = combine_c2f(c2f, tmp,felem_idx,relem_idx_x);
0217                       end
0218                       progress = progress + 1;
0219                       progress_msg(progress,max_iter);
0220                    end
0221                 else
0222                    [fmdl, rmdl, opt, felem_idx] = crop_models(fmdl0,rmdl0,opt,relem_idx_y);
0223                    if any(felem_idx)
0224                       eidors_msg('log_level',eidors_msg('log_level')-2);
0225                       tmp = separable_calculations(fmdl,rmdl,opt);
0226                       eidors_msg('log_level',eidors_msg('log_level')+2);
0227                       c2f = combine_c2f(c2f, tmp,felem_idx,relem_idx_y);
0228                    end
0229                    progress = progress + 1;
0230                    progress_msg(progress, max_iter);
0231                 end
0232              end
0233 
0234           else
0235              [fmdl, rmdl, opt, felem_idx] = crop_models(fmdl0,rmdl0,opt,relem_idx);
0236              if any(felem_idx)
0237                 eidors_msg('log_level',eidors_msg('log_level')-2);
0238                 tmp = separable_calculations(fmdl,rmdl,opt);
0239                 eidors_msg('log_level',eidors_msg('log_level')+2);
0240                 c2f = combine_c2f(c2f, tmp,felem_idx,relem_idx);
0241              end
0242              progress = progress + 1;
0243              progress_msg(progress, max_iter);
0244           end
0245           
0246           
0247        end
0248        progress_msg(Inf);
0249     end
0250 
0251 %-------------------------------------------------------------------------%
0252 % Wrapper helper, combines partial c2f matrices
0253 function c2f = combine_c2f(c2f, tmp,felem_idx,relem_idx)
0254     fidx = find(felem_idx); 
0255     ridx = find(relem_idx);
0256     c2f(fidx,ridx) = c2f(fidx,ridx) + tmp;
0257     
0258 %-------------------------------------------------------------------------%
0259 % The main function
0260 function [c2f, m] = separable_calculations(fmdl,rmdl,opt)
0261     DEBUG = eidors_debug('query','mk_grid_c2f');
0262    
0263     progress_msg('Prepare vox model...');
0264     m = prepare_vox_mdl(rmdl,opt);
0265     progress_msg(Inf);
0266     
0267     try
0268     
0269     % tet edge v. vox face
0270     progress_msg('Find tet_edge2vox_face intersections...')
0271     [intpts1, rec2tedge, rec2intpt1, tedge2intpt1] = get_voxel_intersection_points(fmdl,m.faces,opt);
0272     progress_msg(sprintf('Found %d', size(intpts1,1)), Inf);
0273     
0274     % vox edge v. tet face
0275     progress_msg('Find vox_edge2tet_face intersections...',0,3)
0276     [intpts2, tri2vedge, tri2intpt2, vedge2intpt2] = get_tet_intersection_points(fmdl,m,opt);
0277     progress_msg(sprintf('Found %d', size(intpts2,1)), Inf);
0278     
0279     % Note: Rather than calculating edge2edge intersections, one could
0280     % include them in one or both of the previous tests. However, it is
0281     % then difficult to guarantee that the numerically border-line cases
0282     % get assigned to all the vox and tets concerned
0283     
0284     % vox edge v. tet edge
0285     pmopt.final_msg = 'none';   
0286     progress_msg('Find tet_edge2vox_edge intersections...',-1,pmopt)
0287     [intpts3, tedge2vedge, tedge2intpt3, vedge2intpt3] =find_edge2edge_intersections(fmdl.edges,fmdl.nodes,m.edges,rmdl.nodes, opt.tol_edge2edge); 
0288     progress_msg(sprintf('Found %d',size(intpts3,1)),Inf);
0289     
0290     % tet node in vox
0291     progress_msg('Find tet_nodes in voxels...')
0292     [fnode2vox] = get_nodes_in_voxels(fmdl,rmdl);
0293     progress_msg(sprintf('Found %d', nnz(fnode2vox)), Inf);
0294     
0295     % vox node in tet
0296     progress_msg('Find vox_nodes in tets...');
0297     vnode2tet = get_nodes_in_tets(fmdl,rmdl, opt);
0298     progress_msg(sprintf('Found %d', nnz(vnode2tet)), Inf);
0299     
0300     % vox contained in tet
0301     progress_msg('Find vox contained in tet...')
0302     vox_in_tet = (m.node2vox' * vnode2tet) == 8;
0303     progress_msg(sprintf('Found %d',nnz(vox_in_tet)), Inf);
0304     
0305     % tet contained in vox
0306     progress_msg('Find tets contained in vox...');
0307     tet_in_vox = (double(fmdl.node2elem') * fnode2vox) == 4;
0308     progress_msg(sprintf('Found %d',nnz(tet_in_vox)), Inf);
0309     
0310     
0311     % tets and vox that intersect
0312     progress_msg('Find total vox v. tet intersections...');
0313     vox2intTet =   m.vox2face * (rec2tedge>0) * fmdl.edge2elem ...
0314                  | m.vox2edge * (tri2vedge>0)' * fmdl.elem2face' ...
0315                  | m.vox2edge * tedge2vedge' * fmdl.edge2elem;
0316     % exclude complete inclusion
0317     vox2intTet = vox2intTet & ~vox_in_tet & ~tet_in_vox';
0318     progress_msg(sprintf('Found %d',nnz(vox2intTet)), Inf);
0319     
0320     catch err
0321        if (strcmp(err.identifier,'MATLAB:nomem'))
0322           msg = sprintf('%s', ...
0323              'Matlab ran out of memory. Consider setting save_memory > 0.\n', ...
0324              'See HELP MK_GRID_C2F for details.');
0325           error('EIDORS:mk_grid_c2f:nomem',msg);
0326        else
0327           rethrow(err);
0328        end
0329     end
0330     
0331     progress_msg('Calculate intersection volumes...');
0332     % sparse logical multiplication doesn't exist
0333     vox2intpt1 = logical(m.vox2face*rec2intpt1)'; 
0334     tet2intpt1 = logical(fmdl.edge2elem'*tedge2intpt1)';
0335 
0336     tet2intpt2 = logical(fmdl.elem2face*tri2intpt2)';
0337     vox2intpt2 = logical(m.vox2edge*vedge2intpt2)';
0338 
0339     tet2intpt3 = logical(fmdl.edge2elem'*tedge2intpt3)';
0340     vox2intpt3 = logical(m.vox2edge*vedge2intpt3)';
0341     
0342     vox_todo = find(sum(vox2intTet,2)>0);
0343     C = []; F = []; V = [];
0344     
0345     id = 0; lvox = length(vox_todo);
0346     mint = ceil(lvox/100);
0347     for v = vox_todo'
0348         id = id+1;
0349         if mod(id,mint)==0, progress_msg(id/lvox); end
0350         tet_todo = find(vox2intTet(v,:));
0351         common_intpts1 = bsxfun(@and,vox2intpt1(:,v), tet2intpt1(:,tet_todo));
0352         common_intpts2 = bsxfun(@and,vox2intpt2(:,v), tet2intpt2(:,tet_todo));
0353         common_intpts3 = bsxfun(@and,vox2intpt3(:,v), tet2intpt3(:,tet_todo));
0354         tet_nodes     = bsxfun(@and,fnode2vox(:,v), fmdl.node2elem(:,tet_todo));
0355         vox_nodes     = bsxfun(@and,vnode2tet(:,tet_todo), m.node2vox(:,v));
0356 
0357         C = [C; v*ones(numel(tet_todo),1)];
0358         F = [F; tet_todo'];
0359         last_v = numel(V);
0360         V = [V; zeros(numel(tet_todo),1)]; % pre-allocate
0361 
0362         for t = 1:numel(tet_todo)
0363             pts = [ intpts1(common_intpts1(:,t),:); % tet edge v. vox face
0364                     intpts2(common_intpts2(:,t),:); % vox edge v. tet face
0365                     intpts3(common_intpts3(:,t),:); % vox edge v. tet edge
0366                     fmdl.nodes(tet_nodes(:,t),:);   % tet node in vox
0367                     rmdl.nodes(vox_nodes(:,t),:)];  % vox node in tet
0368             last_v = last_v + 1;
0369             ok = false;
0370             if size(pts,1) < 4 
0371               % we should have enough confidence by now
0372               % unless the tolerances were specified weird
0373               continue 
0374             end
0375             if size(pts,1) < 4 % test if edge lies on the plane of the vox
0376                 % check for edges along the x y or z axis
0377                 % this includes coplanar faces
0378                 E = fmdl.edges(fmdl.elem2edge(tet_todo(t),:),:);
0379                 P1 = fmdl.nodes(E(:,1),:);
0380                 P2 = fmdl.nodes(E(:,2),:);
0381                 % this test is sensitive, but not specific
0382                 % it should also check if both pts come from the same edge and
0383                 % that edge fullfils the condition
0384                 D = P1-P2;
0385                 ok = any(abs(D(:)) <= eps); 
0386             end 
0387             if ok; continue, end % otherwise convhulln will throw an error
0388             try
0389                 % move points to origin (helps for small elements at
0390                 % large coordinates
0391                 ctr = mean(pts);
0392                 pts = bsxfun(@minus,pts,ctr);
0393                 scale = max(abs(pts(:)));
0394                 if scale == 0 %happens when there's only one point
0395                    continue
0396                 end
0397                 % scale largest coordinate to 1 (helps with precision)
0398                 pts = pts ./ scale;
0399                 % force thorough search for initinal simplex and
0400                 % supress precision warnings
0401                 [~, V(last_v)] = convhulln_clean(pts);
0402                 V(last_v) = V(last_v) * scale^3; % undo scaling
0403             catch err
0404                 ok = false;
0405                 switch err.identifier
0406                     case {'MATLAB:qhullmx:DegenerateData'}
0407                         % if QHull is degenerate, then are is zero.
0408                         ok = 1;
0409 
0410                     case {'MATLAB:qhullmx:DegenerateData', 'MATLAB:qhullmx:UndefinedError', ...
0411                             'MATLAB:cgprechecks:NotEnoughPts'}
0412                         % check for edges along the x y or z axis
0413                         % this includes coplanar faces
0414                         E = fmdl.edges(fmdl.elem2edge(tet_todo(t),:),:);
0415                         P1 = fmdl.nodes(E(:,1),:);
0416                         P2 = fmdl.nodes(E(:,2),:);
0417                         % this test is sensitive, but not specific
0418                         % it should also check if both pts come from the same edge and
0419                         % that edge fullfils the condition
0420                         D = P1-P2;
0421                         ok = any(abs(D) < 10*eps);
0422                         % edge-edge intersections often appear to also
0423                         % cross faces, there doesn't seem to be a good
0424                         % specific way to catch that
0425                         u = uniquetol(pts*scale,10*eps,'ByRows',true,'DataScale', 1);
0426                         ok = ok | size(u,1) < 4;
0427                 end
0428                 if ~ok & size(u,1) == 4
0429                         ok = ok | det([ones(4,1),u])<eps;
0430                 end
0431                 if ~ok
0432                     if DEBUG || eidors_debug('query','mk_grid_c2f:convhulln')
0433                         tet.nodes = fmdl.nodes;
0434                         vox.nodes = rmdl.nodes;
0435                         tet.type = 'fwd_model';
0436                         vox.type = 'fwd_model';
0437                         vox.elems = m.faces(logical(m.vox2face(v,:)),:);
0438                         vox.boundary = vox.elems;
0439                         tet.elems = fmdl.elems(tet_todo(t),:);
0440                         
0441                         clf
0442                         show_fem(vox)
0443                         hold on
0444                         h = show_fem(tet);
0445                         set(h,'EdgeColor','b')
0446                         pts = bsxfun(@plus,pts*scale,ctr);
0447                         plot3(pts(:,1),pts(:,2),pts(:,3),'o');
0448 %                         plot3(nt(:,1),nt(:,2),nt(:,3),'xr');
0449 %                         plot3(nv(:,1),nv(:,2),nv(:,3),'xb');
0450                         hold off
0451                         axis auto
0452                         keyboard
0453                     else
0454                         fprintf('\n');
0455                         eidors_msg(['convhulln has thrown an error. ' ...
0456                             'Enable eidors_debug on mk_grid_c2f:convhulln and re-run to see a debug plot'],0);
0457                         rethrow(err);
0458                     end
0459                 end
0460             end
0461         end
0462     end
0463     progress_msg(Inf);
0464     c2f = sparse(F,C,V,opt.nTet,opt.nVox);
0465     
0466     % add vox contained in tet
0467     c2f = c2f + bsxfun(@times, sparse(vox_in_tet), m.volume)';
0468     
0469     % normalize to tet volume
0470     vol = get_elem_volume(fmdl, 'no_c2f');
0471     c2f = bsxfun(@rdivide,c2f,vol);
0472 
0473     % add tets contained in vox
0474 
0475     c2f = c2f + tet_in_vox;
0476 
0477 %-------------------------------------------------------------------------%
0478 % Crop obviously non-overlapping parts of the models to limit computation
0479 function [fmdl, rmdl, opt, felem_idx] = crop_models(fmdl0,rmdl0,opt, relem_idx)
0480    
0481    fmdl = fmdl0;
0482    rmdl = rmdl0;
0483    
0484    opt.nVox = opt.Xsz*opt.Ysz*opt.Zsz;
0485    idx = rmdl0.coarse2fine * relem_idx > 0;
0486    rmdl.elems = rmdl0.elems(idx,:);
0487    [node_idx, m,n] = unique(rmdl.elems(:));
0488    rmdl.elems = reshape(n,size(rmdl.elems));
0489    rmdl.nodes = rmdl0.nodes(node_idx,:);
0490    
0491    fnode_above = fmdl0.nodes(:,3) > opt.zvec(end);
0492    % matlab is so stupid!
0493    felem_above = all(reshape(fnode_above(fmdl0.elems), size(fmdl0.elems)),2);
0494    
0495    fnode_below = fmdl0.nodes(:,3) < opt.zvec(1);
0496    felem_below = all(reshape(fnode_below(fmdl0.elems), size(fmdl0.elems)),2);
0497    
0498    felem_idx   = ~(felem_below | felem_above);
0499    fmdl.elems = fmdl0.elems(felem_idx,:);
0500    fnode_idx = unique(fmdl.elems(:));
0501    fnode_idx_map = zeros(size(fmdl0.nodes,1),1);
0502    fnode_idx_map(fnode_idx) = 1:length(fnode_idx);
0503    fmdl.elems = reshape(fnode_idx_map(fmdl.elems),size(fmdl.elems));
0504    fmdl.edges = reshape(fnode_idx_map(fmdl0.edges),size(fmdl0.edges));
0505    fmdl.nodes = fmdl0.nodes(fnode_idx,:);
0506    fmdl.node2elem = fmdl0.node2elem(fnode_idx,felem_idx);
0507    
0508    
0509    fface_idx = sum(fmdl0.elem2face(felem_idx,:),1)>0;
0510    fmdl.elem2face = fmdl0.elem2face(felem_idx,fface_idx);
0511    felem_idx_map = zeros(size(felem_idx));
0512    felem_idx_map(felem_idx) = 1:nnz(felem_idx);
0513    felem_idx_map = [0; felem_idx_map];
0514    fmdl.face2elem = fmdl0.face2elem(fface_idx,:);
0515    fmdl.face2elem = reshape(felem_idx_map(fmdl.face2elem + 1),size(fmdl.face2elem));
0516    fmdl.normals = fmdl0.normals(fface_idx,:);
0517    fmdl.faces = fmdl0.faces(fface_idx,:);
0518    fmdl.faces = reshape(fnode_idx_map(fmdl.faces),size(fmdl.faces));
0519    
0520    fedge_idx = unique(fmdl0.elem2edge(felem_idx,:));
0521    fedge_idx_map = zeros(size(fmdl0.edges,1),1);
0522    fedge_idx_map(fedge_idx) = 1:length(fedge_idx);
0523    fmdl.elem2edge = fedge_idx_map(fmdl0.elem2edge(felem_idx,:));
0524    fmdl.edge2elem = fmdl0.edge2elem(fedge_idx,felem_idx);
0525    fmdl.edges = fmdl.edges(fedge_idx,:);
0526    
0527    opt.nTet = size(fmdl.elems,1);
0528 
0529 %-------------------------------------------------------------------------%
0530 % Prepare matrices for the voxel model
0531 function fmdl = prepare_fmdl(fmdl)
0532     fmopt.elem2edge = true;
0533     fmopt.edge2elem = true;
0534     fmopt.face2elem = true;
0535     fmopt.node2elem = true;
0536     fmopt.normals   = true;
0537     fmopt.linear_reorder = false; % this is slow and not needed
0538     ll = eidors_msg('log_level',1);
0539     fmdl = fix_model(fmdl,fmopt);
0540     eidors_msg('log_level',ll);
0541     fmdl.node2elem = logical(fmdl.node2elem);
0542     nElem = size(fmdl.elems,1);
0543     nFace = size(fmdl.faces,1);
0544     fmdl.elem2face = sparse(repmat((1:nElem)',1,4),double(fmdl.elem2face),true,nElem,nFace);
0545 
0546 
0547 %-------------------------------------------------------------------------%
0548 % Prepare matrices for the voxel model
0549 function m = prepare_vox_mdl(rmdl,opt)
0550 
0551     DEBUG = eidors_debug('query','mk_grid_c2f');
0552 
0553     [voxels, m.node2vox] = mk_voxels(opt);
0554     if DEBUG || eidors_debug('query','mk_grid_c2f:mk_voxels')
0555         show_voxels(rmdl,voxels); title('mk\_voxels');
0556     end
0557 
0558     m.faces = mk_faces(voxels,opt);
0559     if DEBUG || eidors_debug('query','mk_grid_c2f:mk_faces')
0560         test_faces(rmdl,m.faces,opt); title('mk\_faces');
0561     end
0562 
0563     m.vox2face = mk_vox2face(opt);
0564     if DEBUG || eidors_debug('query','mk_grid_c2f:mk_vox2face')
0565         % the numbers shown are useless, just check if all faces are present
0566         show_voxels(rmdl,m.faces(any(m.vox2face),:));title('mk\_vox2face');
0567     end
0568     m.edges = mk_edges(voxels,opt);
0569 
0570     m.edge_length = rmdl.nodes(m.edges(:,1),:) - rmdl.nodes(m.edges(:,2),:);
0571     m.edge_length = sqrt(sum(m.edge_length.^2,2));
0572     
0573     [m.vox2edge, m.volume]= mk_vox2edge(m,opt);
0574 
0575 
0576 %-------------------------------------------------------------------------%
0577 % Assign each rmdl node to the tet it is in (nodes on tet faces are counted
0578 % mutltiple times)
0579 function rnode2tet = get_nodes_in_tets(fmdl,rmdl, opt)
0580     
0581     [A,b] = tet_to_inequal(fmdl.nodes,fmdl.elems);
0582     progress_msg(.01);
0583     % This is split to decrease the memory footprint
0584     rnode2tet = (bsxfun(@minus, A(1:4:end,:)*rmdl.nodes',b(1:4:end)) <= opt.tol_node2tet)';
0585     progress_msg(.21);
0586     for i = 2:4 
0587         rnode2tet = rnode2tet & (bsxfun(@minus, A(i:4:end,:)*rmdl.nodes',b(i:4:end)) <= opt.tol_node2tet)'; 
0588         progress_msg(.21 + (i-1)*.23);
0589     end
0590     
0591     % exclude coinciding nodes
0592     ex= bsxfun(@eq,rmdl.nodes(:,1),fmdl.nodes(:,1)') & ...
0593         bsxfun(@eq,rmdl.nodes(:,2),fmdl.nodes(:,2)') & ...
0594         bsxfun(@eq,rmdl.nodes(:,3),fmdl.nodes(:,3)');
0595     progress_msg(.94);
0596     rnode2tet(any(ex,2),:) = 0;
0597     progress_msg(1);
0598     
0599 
0600 %-------------------------------------------------------------------------%
0601 % Assign each fmdl node to the vox it is in (nodes on vox faces are counted
0602 % mutltiple times)
0603 function [insnode] = get_nodes_in_voxels(fmdl,rmdl)
0604 
0605     E = reshape(rmdl.elems',4*6,[])';
0606     E = E(:,[1 2 3 4 8 12 16 23]);
0607 
0608     NE = size(E,1);
0609     xnodes = reshape(rmdl.nodes(E,1),NE,[]);
0610     ynodes = reshape(rmdl.nodes(E,2),NE,[]);
0611     znodes = reshape(rmdl.nodes(E,3),NE,[]);
0612     minx = min(xnodes,[],2);
0613     maxx = max(xnodes,[],2);
0614     miny = min(ynodes,[],2);
0615     maxy = max(ynodes,[],2);
0616     minz = min(znodes,[],2);
0617     maxz = max(znodes,[],2);
0618 
0619     leftof  = bsxfun(@lt, fmdl.nodes(:,1)+eps, minx');
0620     rightof = bsxfun(@gt, fmdl.nodes(:,1)-eps, maxx');
0621     infront = bsxfun(@lt, fmdl.nodes(:,2)+eps, miny');
0622     behind  = bsxfun(@gt, fmdl.nodes(:,2)-eps, maxy');
0623     below   = bsxfun(@lt, fmdl.nodes(:,3)+eps, minz');
0624     above   = bsxfun(@gt, fmdl.nodes(:,3)-eps, maxz');
0625 
0626     outnode = leftof | rightof | behind | infront | below | above;
0627     insnode = sparse(~outnode);
0628 
0629 %-------------------------------------------------------------------------%
0630 % Calculate intersection points between vox faces and tet edges
0631 function [intpts, face2edge, face2intpt, edge2intpt] = get_voxel_intersection_points(fmdl,faces,opt)
0632 edges = fmdl.edges;
0633 nodes = fmdl.nodes;
0634 dir   = nodes(edges(:,2),:) - nodes(edges(:,1),:);
0635 intpts = [];
0636 F = []; E = []; I = [];
0637 
0638 SZ = [opt.Xsz, opt.Ysz, opt.Zsz];
0639 VEC = {opt.xvec, opt.yvec, opt.zvec};
0640 % show_voxels(rmdl, voxels(faces,:))
0641 
0642 %  mdl = rmdl;
0643 %  mdl.elems = voxels(faces,:);
0644 %  img = mk_image(mdl,1);
0645 todo = sum(SZ)+3;
0646 id = 0;
0647 step = ceil(todo/100);
0648 
0649 for d = 1:3
0650     for x = 0:SZ(d)
0651         id = id+1;
0652         if mod(id,step)==0, progress_msg(id/todo); end
0653         % plane-edge intersection
0654         t = (VEC{d}(x+1) - nodes(edges(:,1),d)) ./ dir(:,d);
0655         crossing = t>0 & t<1;
0656         crossed  = find(crossing);
0657         if isempty(crossed), continue, end;
0658         % intersection points
0659         tmp   = nodes(edges(crossing,1),:) + bsxfun(@times,t(crossing),dir(crossing,:));
0660         face_coord = zeros(length(crossed),3);
0661         rmv = false(length(crossed),1);
0662         for i = 1:3
0663             if i == d
0664                face_coord(:,i) = x+1;
0665             else
0666                 face_coord(:,i) = sum(bsxfun(@times,diff(bsxfun(@lt, tmp(:,i), VEC{i}'),1,2),(1:SZ(i))),2);
0667                 rmv = rmv | face_coord(:,i) == 0;
0668             end
0669         end
0670         % also reject intersections with voxel nodes and edges
0671         same_x = any(bsxfun(@eq,tmp(:,1),opt.xvec'),2);
0672         same_y = any(bsxfun(@eq,tmp(:,2),opt.yvec'),2);
0673         same_z = any(bsxfun(@eq,tmp(:,3),opt.zvec'),2);
0674         rmv = rmv | (same_x + same_y + same_z) > 1;
0675         if nnz(rmv) == numel(rmv), continue, end
0676         tmp(rmv,:) = []; 
0677         face_coord(rmv,:) = [];
0678         crossed(rmv) = [];
0679         face_coord = face_coord - 1;
0680         face_idx = d + face_coord(:,1)*opt.xstep + face_coord(:,2)*opt.ystep + face_coord(:,3)*opt.zstep;
0681         I = [I; (1:size(tmp,1))' + size(I,1)];
0682         F = [F; face_idx];
0683         E = [E; crossed];
0684         intpts = [intpts; tmp];
0685 %         unique(face_idx)
0686 %         img.elem_data(:) = 0;
0687 %         img.elem_data(unique(face_idx)) = 1;
0688 %         show_fem(img,[0 0 1]);
0689         %     hold on
0690         %     x1 =fmdl.nodes(fmdl.edges(crossing,1),1);
0691         %     x2 =fmdl.nodes(fmdl.edges(crossing,2),1);
0692         %     y1 =fmdl.nodes(fmdl.edges(crossing,1),2);
0693         %     y2 =fmdl.nodes(fmdl.edges(crossing,2),2);
0694         %     z1 =fmdl.nodes(fmdl.edges(crossing,1),3);
0695         %     z2 =fmdl.nodes(fmdl.edges(crossing,2),3);
0696         %     plot3([x1 x2]',[y1 y2]',[z1 z2]','b');
0697         %     plot3(intpts(:,1),intpts(:,2),intpts(:,3),'ro','MarkerFaceColor','r');
0698         %     hold off
0699         %     axis auto
0700         %     view(2)
0701         %     pause
0702     end
0703 end
0704 face2edge = sparse(F,E,I,size(faces,1),size(edges,1));
0705 face2intpt = sparse(F,I,ones(size(I)),size(faces,1),size(I,1));
0706 edge2intpt  = sparse(E,I,ones(size(I)),size(edges,1),size(I,1));
0707 
0708 %-------------------------------------------------------------------------%
0709 % Calculate intersection points between tet faces and vox edges
0710 function [intpts, tri2edge, tri2intpt, edge2intpt] = get_tet_intersection_points(fmdl,m,opt)
0711     
0712     intpts = [];
0713     T = []; E = []; I = [];
0714     Xsz = opt.Xsz; Ysz = opt.Ysz; Zsz = opt.Zsz;
0715     SZ = [opt.Xsz, opt.Ysz, opt.Zsz];
0716     VEC = {opt.xvec, opt.yvec, opt.zvec};
0717     STEP(1) = opt.xstep; 
0718     STEP(2) = STEP(1)*(Xsz+1);
0719     STEP(3) = STEP(2)*(Ysz+1);
0720     
0721     d = sum(fmdl.normals .* fmdl.nodes(fmdl.faces(:,1),:),2);
0722 
0723     line_axis = [3 1 2];
0724     for i = 1:3
0725         progress_msg(i,3);
0726         % define edge lines
0727         idx = 1:3;
0728         op = line_axis(i);
0729         idx(op) = [];
0730         
0731         pts = [repmat(VEC{idx(1)},SZ(idx(2))+1,1) kron(VEC{idx(2)},ones(SZ(idx(1))+1,1)) ];
0732         pt_idx = uint32(repmat((0:SZ(idx(1)))',SZ(idx(2))+1,1)*STEP(idx(1)) ...
0733                 + kron((0:SZ(idx(2)))', ones(SZ(idx(1))+1,1))*STEP(idx(2)));
0734         
0735         % project on faces
0736         % plane equation is ax+by+cz+d = 0, where d = -(ax0 + by0 + cz0)
0737         z = repmat(d,1,size(pts,1));
0738         for j = 1:2
0739             z = z - bsxfun(@times,fmdl.normals(:,idx(j)),pts(:,j)');
0740         end
0741         z = z ./ repmat(fmdl.normals(:,op),1,size(pts,1));
0742         in = point_in_triangle(pts,fmdl.faces,fmdl.nodes(:,idx),2*opt.tol_edge2tri)';
0743 
0744         
0745         % reject voxel nodes
0746         v = VEC{op}';
0747         in = in & ~reshape(any(bsxfun(@eq,z(:),v),2),size(in));
0748         M = bsxfun(@lt, z(:),v);
0749         M = xor(M(:,1:end-1), M(:,2:end));
0750 %       edge_num1= reshape(uint32(sum(bsxfun(@times,uint32(M),uint32(1:SZ(op))),2)), size(z));
0751         edge_num = reshape(uint32(M*(1:SZ(op))'), size(z));
0752 %       if ~(all(all(edge_num1==edge_num))); keyboard; end
0753         in = in & edge_num > 0;
0754         if nnz(in) == 0, continue, end
0755         edge_num(~in) = 0;
0756 
0757         edge_num(in) = (edge_num(in)-1) * STEP(op);
0758         edge_idx = edge_num + bsxfun(@times,uint32(full(in)), uint32(i) + pt_idx');
0759 
0760         [t, p] = find(in);
0761         tmp = zeros(length(p),3);
0762         tmp(:,idx) = pts(p,1:2);
0763         tmp(:,op) = z(in);
0764         
0765         I = [I; (1:size(tmp,1))' + size(I,1)];
0766         T = [T; t];
0767         E = [E; edge_idx(in)];
0768         intpts = [intpts; tmp];        
0769     end
0770     
0771     E = double(E); % keep sparse happy
0772     tri2edge = sparse(T,E,I,size(fmdl.faces,1),size(m.edges,1));
0773     tri2intpt = sparse(T,I,ones(size(I)),size(fmdl.faces,1),size(I,1));
0774     edge2intpt  = sparse(E,I,ones(size(I)),size(m.edges,1),size(I,1));
0775 
0776    
0777 %-------------------------------------------------------------------------%
0778 % Make voxels
0779 function [voxels, node2vox] = mk_voxels(opt)
0780     Xsz = opt.Xsz; Ysz = opt.Ysz; Zsz = opt.Zsz;
0781     Xp = Xsz+1; Yp = Ysz+1; Zp = Zsz+1; % number of planes
0782 
0783     up = Xp*Yp;
0784     vox = [ 1       1+up    1+up+Xp     1+Xp;
0785             1       2       2+up        1+up;
0786             1       1+Xp    2+Xp        2;
0787             2       2+up    2+up+Xp     2+Xp;
0788             1+Xp    2+Xp    2+up+Xp     1+up+Xp;
0789             1+up    1+Xp+up 2+Xp+up     2+up];
0790 
0791     voxrow   = bsxfun(@plus, repmat(vox,     Xsz,1), kron(0:Xsz-1          ,ones(1,6))');
0792     voxplane = bsxfun(@plus, repmat(voxrow,  Ysz,1), kron(Xp*(0:Ysz-1) , ones(1, 6*Xsz))');
0793     voxels   = bsxfun(@plus, repmat(voxplane,Zsz,1), kron(up*(0:Zsz-1) , ones(1, 6*Xsz*Ysz))');
0794     nVox     = Xsz * Ysz * Zsz;
0795     node2vox = sparse(voxels(1:3:end,:),kron((1:nVox)',ones(2,4)),1);
0796     
0797 %-------------------------------------------------------------------------%
0798 % Make faces
0799 function faces = mk_faces(voxels,opt)
0800     Xsz = opt.Xsz; Ysz = opt.Ysz; Zsz = opt.Zsz;
0801     
0802     facerow     = [nonzeros(bsxfun(@plus,(1:3)',0:6:(Xsz-1)*6)); (Xsz-1)*6+4];
0803     faceplane   = bsxfun(@plus,facerow,0:6*Xsz:6*Xsz*(Ysz-1));
0804     % cap duplicates faces in order to preserve nice indexing
0805     cap         = kron(ones(3,1),faceplane(2:3:end,end)')+3;
0806     faceplane   = [faceplane,[ cap(:); cap(end)]];
0807     faces       = bsxfun(@plus, faceplane(:), 0:6*Xsz*Ysz:6*Xsz*Ysz*(Zsz-1));
0808     % cap duplicates faces in order to preserve nice indexing
0809     cap         = reshape(kron(ones(3,1), 3:3:3*Xsz*Ysz'),3*Xsz,[]);
0810     cap         = bsxfun(@plus, cap, 0:Ysz-1);
0811     cap(end+1,:)= cap(end,:);
0812     faces       = [faces(:); faces(cap(:),end)+3];
0813     faces       = voxels(faces,:);
0814     
0815 %-------------------------------------------------------------------------%
0816 % Make edges
0817 function edges = mk_edges(voxels, opt)
0818     Xsz = opt.Xsz; Ysz = opt.Ysz; Zsz = opt.Zsz;
0819     
0820     edgerow     = voxels([nonzeros(bsxfun(@plus,(1:3)',0:6:(Xsz-1)*6)); (Xsz-1)*6+4],1:2);
0821     edgerow(end+1,:) = edgerow(end,:);
0822     edgerow(end+1,:) = voxels((Xsz-1)*6+4, [1 4]);
0823     edgeplane = repmat(edgerow,Ysz+1,1) + kron((Xsz+1)*(0:Ysz)', ones(size(edgerow)));
0824     % replace ficticious edges with repetitions;
0825     edgeplane((size(edgerow,1)*Ysz +3):3:end,:) = edgeplane((size(edgerow,1)*Ysz +2):3:end,:);
0826 
0827     up = (Xsz+1)*(Ysz+1);
0828     edges = repmat(edgeplane,Zsz+1,1) + kron((0:Zsz)'*up, ones(size(edgeplane)));
0829     
0830     % replace ficticious edges with repetitions;
0831     start = (size(edgeplane,1)*Zsz);
0832     edges( (start +1):3:end,:) = edges( (start +2):3:end,:);
0833     start = (size(edgeplane,1)*Zsz) + 3*Xsz;
0834     step  = size(edgerow,1);
0835     edges( (start +1):step:end,:) = edges( (start +3):step:end,:);
0836     edges( (start +2):step:end,:) = edges( (start +3):step:end,:);
0837     edges( end-2:end,:) = edges(end-5:end-3,:);
0838 
0839 %-------------------------------------------------------------------------%
0840 % Make mapping between voxels and faces
0841 function vox2face = mk_vox2face(opt)
0842     Xsz = opt.Xsz; Ysz = opt.Ysz; Zsz = opt.Zsz;
0843     xstep = opt.xstep; ystep = opt.ystep; zstep = opt.zstep;
0844     
0845     vox = [1 2 3 1+xstep 2+ystep 3+zstep];
0846     voxrow = bsxfun(@plus,(0:Xsz-1)'*xstep,vox)';
0847     voxplane = bsxfun(@plus,(0:Ysz-1)*ystep,voxrow(:));
0848     voxvol   = bsxfun(@plus,(0:Zsz-1)*zstep,voxplane(:));
0849     voxvol   = reshape(voxvol, 6, [])';
0850     I = kron((1:size(voxvol,1))',ones(1,6));
0851     nFaces  = Zsz*(Ysz+1)*(3*Xsz+1) + Ysz*(3*Xsz+1);
0852     nVox    = Xsz*Ysz*Zsz;
0853     vox2face = sparse(I,voxvol,ones(size(I)),nVox,nFaces);   
0854     
0855 %-------------------------------------------------------------------------%
0856 % Make mapping between voxels and edges
0857 function [vox2edge, vol] = mk_vox2edge(m,opt)
0858     Xsz = opt.Xsz; Ysz = opt.Ysz; Zsz = opt.Zsz;
0859     xstep = opt.xstep; ystep = xstep*(Xsz+1); zstep = ystep*(Ysz+1);
0860     
0861     vox = [1 2 3 1+xstep 3+xstep 1+ystep 2+ystep 1+ystep+xstep ...
0862            2+zstep 3+zstep 3+zstep+xstep 2+zstep+ystep];
0863     voxrow = bsxfun(@plus,(0:Xsz-1)'*xstep,vox)';
0864     voxplane = bsxfun(@plus,(0:Ysz-1)*ystep,voxrow(:));
0865     voxvol   = bsxfun(@plus,(0:Zsz-1)*zstep,voxplane(:));
0866     voxvol   = reshape(voxvol, 12, [])';
0867     I = kron((1:size(voxvol,1))',ones(1,12));
0868     nEdges  = 3*(Zsz+1)*(Ysz+1)*(Xsz+1);
0869     nVox    = Xsz*Ysz*Zsz;
0870     vox2edge = sparse(I,voxvol,ones(size(I)),nVox,nEdges);
0871     vol =      m.edge_length(voxvol(:,1)) ...
0872             .* m.edge_length(voxvol(:,2)) ...
0873             .* m.edge_length(voxvol(:,3));
0874     
0875     
0876 %-------------------------------------------------------------------------%
0877 % Show voxels
0878 function show_voxels(rmdl,voxels)    
0879     mdl=rmdl;
0880     mdl.elems = voxels;
0881     mdl.boundary = mdl.elems;
0882     clf
0883     show_fem(mdl,[0 0 1]);
0884 
0885 %-------------------------------------------------------------------------%
0886 % Test indexing of faces on planes by showing random x, y and z plane
0887 function test_faces(rmdl, faces, opt)
0888     mdl = rmdl;
0889     mdl.elems = faces;
0890     mdl.boundary = mdl.elems;
0891     img = mk_image(mdl,0);
0892     img.elem_data(opt.xplane + (randi(opt.Xsz+1)-1)*opt.xstep) = 1;
0893     img.elem_data(opt.yplane + (randi(opt.Ysz+1)-1)*opt.ystep) = 2;
0894     img.elem_data(opt.zplane + (randi(opt.Zsz+1)-1)*opt.zstep) = 3;
0895     show_fem(img,[0 0 1]);
0896 
0897 %-------------------------------------------------------------------------%
0898 % Center scale models
0899 function [fmdl,rmdl] = center_scale_models(fmdl,rmdl, opt)
0900     % ensure 3D
0901     if size(rmdl.nodes,2)==2;
0902        rmdl = rmdl_3D(rmdl, fmdl);
0903     end
0904     ctr = mean([min(rmdl.nodes);max(rmdl.nodes)]);
0905     rmdl.nodes = bsxfun(@minus,rmdl.nodes,ctr);
0906     fmdl.nodes = bsxfun(@minus,fmdl.nodes,ctr);
0907     if isfield(opt,'do_not_scale') && opt.do_not_scale    
0908         return
0909     end
0910     maxnode = min( max(abs(rmdl.nodes(:))), max(abs(fmdl.nodes(:))));
0911     scale = 1/maxnode;
0912     rmdl.nodes = scale*rmdl.nodes;
0913     fmdl.nodes = scale*fmdl.nodes;
0914     eidors_msg('@@ models scaled by %g', scale,2);
0915     
0916 
0917 % Make a 3D grid from 2D model
0918 % where the max and min enclose fmdl
0919 function rmdl = rmdl_3D(rmdl, fmdl)
0920    fzmax = max(fmdl.nodes(:,3));
0921    fzmin = min(fmdl.nodes(:,3));
0922    rm_nodes = rmdl.nodes;
0923    oo = 0*rm_nodes(:,1) + 1;
0924    N = num_nodes(rmdl);
0925    % nodes at top + bottom
0926    rmdl.nodes = [ ...
0927       rmdl.nodes,oo*fzmin;
0928       rmdl.nodes,oo*fzmax];
0929    rme = [
0930       rmdl.elems(:,[1,2,3,1]) + [N,0,0,0];
0931       rmdl.elems(:,[1,2,3,2]) + [N,N,0,0];
0932       rmdl.elems(:,[1,2,3,3]) + [N,N,N,0]];
0933    % put elem block together
0934    rme = reshape(rme,[],3,4);
0935    rme = permute(rme,[2,1,3]);
0936    rmdl.elems = reshape(rme,[],4);
0937    rmdl.coarse2fine = kron( ...
0938         rmdl.coarse2fine,[1;1;1]);
0939    return
0940 % Test code to see block
0941    x = ones(6,1);
0942    rmdl.boundary = find_boundary(rmdl);
0943    x(3) = 2;
0944    img = mk_image(rmdl,x);
0945    show_fem(img);
0946        
0947 %-------------------------------------------------------------------------%
0948 % Parse option struct
0949  function opt = parse_opts(fmdl,rmdl, opt)
0950     opt.xvec = unique(rmdl.nodes(:,1));
0951     opt.yvec = unique(rmdl.nodes(:,2));
0952     opt.zvec = unique(rmdl.nodes(:,3));
0953     opt.Xsz  = numel(opt.xvec)-1;
0954     opt.Ysz  = numel(opt.yvec)-1;
0955     opt.Zsz  = numel(opt.zvec)-1;
0956     opt.xstep = 3;
0957     opt.ystep = opt.xstep*opt.Xsz+1;
0958     opt.zstep = opt.ystep*(opt.Ysz+1);
0959     xrow    = 1 + (0:opt.Ysz-1)*opt.ystep;
0960     opt.xplane  = bsxfun(@plus, (0:opt.Zsz-1)'*opt.zstep,xrow);
0961     yrow    = 2 + (0:opt.Xsz-1)*opt.xstep;
0962     opt.yplane  = bsxfun(@plus, (0:opt.Zsz-1)'*opt.zstep,yrow);
0963     zrow    = 3 + (0:opt.Xsz-1)*opt.xstep;
0964     opt.zplane  = bsxfun(@plus, (0:opt.Ysz-1)'*opt.ystep,zrow);
0965     opt.nVox = opt.Xsz*opt.Ysz*opt.Zsz;
0966     opt.nTet = size(fmdl.elems,1);
0967     
0968     if ~isfield(opt, 'tol_node2tet');
0969         opt.tol_node2tet = eps; % * max(rmdl_rng,fmdl_rng)^3;
0970     end
0971     if ~isfield(opt, 'tol_edge2edge')
0972         opt.tol_edge2edge = 6*sqrt(3)*eps(min(max(abs(fmdl.nodes(:))),max(abs(rmdl.nodes(:)))));
0973     end
0974     if ~isfield(opt, 'tol_edge2tri')
0975         opt.tol_edge2tri = eps; %1e-10
0976     end
0977     if ~isfield(opt, 'save_memory')
0978        opt.save_memory = 0;
0979     end
0980     eidors_msg('@@ node2tet  tolerance = %g', opt.tol_node2tet,3);
0981     eidors_msg('@@ edge2edge tolerance = %g', opt.tol_edge2edge,3);
0982     eidors_msg('@@ edge2tri  tolerance = %g', opt.tol_edge2tri,3);
0983 
0984 %-------------------------------------------------------------------------%
0985 % fprintf wrapper to use eidors log_level
0986 function logmsg(varargin)
0987     if eidors_msg('log_level') >= 2
0988         fprintf(varargin{:});
0989     end
0990     
0991    
0992     
0993 %-------------------------------------------------------------------------%
0994 % Perfom unit tests
0995 function do_unit_test
0996     mk_grid_c2f('save_memory',0); do_small_test;
0997     mk_grid_c2f('save_memory',1); do_small_test;
0998     mk_grid_c2f('save_memory',10); do_small_test;
0999     clf
1000     do_case_tests;
1001     return % the tests below are more about performance than correctness
1002     do_realistic_test;
1003     do_edge2edge_timing_test;
1004     
1005 %-------------------------------------------------------------------------%
1006 % Test individual intersections
1007 function do_case_tests
1008     ll = eidors_msg('log_level');
1009     eidors_msg('log_level',1);
1010     vox = mk_grid_model([],0:1,0:1,0:1);
1011     tet.type = 'fwd_model';
1012     tet.elems = [1 2 3 4];
1013 
1014     X = 4; Y = 6;
1015     for i = 1:30
1016         tet.nodes = [0 0 0; 0 1 0; 1 0 0; 0 0 1];
1017         fprintf('%d\n',i);
1018         switch i
1019             case 1 % nothing in common
1020                 txt = 'Nothing in common';
1021                 tet.nodes = tet.nodes + 2;
1022                 subplot(X,Y,i), show_test(vox,tet);
1023                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1024                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1025                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1026                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1027                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1028                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1029                 mk_grid_c2f(tet,vox);
1030 
1031             case 2 % common node
1032                 tet.nodes = tet.nodes + 1;
1033                 subplot(X,Y,i), show_test(vox,tet);
1034                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1035                 unit_test_cmp('Comm node tedge2rec',size(intpts,1),0,0); % nothing
1036                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1037                 unit_test_cmp('Comm node vedge2tri',size(intpts,1),0,0); % nothing
1038                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1039                 unit_test_cmp('Comm node edge2edge',size(intpts,1),0,0); % nothing
1040                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1041                 unit_test_cmp('Comm node tnode2vox',fnode2vox,[1;0;0;0],0); % 1 point
1042                 rnode2tet = test_vnode_in_tet(vox,tet);
1043                 unit_test_cmp('Comm node vnode2tet',nnz(rnode2tet),0,0);  % nothing
1044                 mk_grid_c2f(tet,vox);
1045 
1046             case 3 % tet_edge v vox_node
1047                 txt = 'tet_edge v vox_node';
1048                 tet.nodes(:,1:2) = tet.nodes(:,1:2) - 0.5;
1049                 subplot(X,Y,i), show_test(vox,tet);
1050                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1051                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1052                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1053                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1054                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1055                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1056                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1057                 unit_test_cmp(txt,nnz(fnode2vox),0,0); % nothing
1058                 rnode2tet = test_vnode_in_tet(vox,tet);
1059                 res = zeros(8,1); res(1) = 1;
1060                 unit_test_cmp(txt,rnode2tet,res,0);  % 1 point
1061                 mk_grid_c2f(tet,vox);
1062 
1063             case 4 % tet_edge v vox_edge
1064                 txt = 'tet_edge v vox_edge';
1065                 tet.nodes(:,1:2) = tet.nodes(:,1:2) - 0.5;
1066                 tet.nodes(:,3)   = tet.nodes(:,3)   + 0.5;
1067                 subplot(X,Y,i), show_test(vox,tet);
1068                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1069                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1070                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1071                 unit_test_cmp(txt,size(unique(intpts,'rows'),1),1,0); % 1 point
1072                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1073                 unit_test_cmp(txt,size(intpts,1),1,0); % 1 point
1074                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1075                 unit_test_cmp(txt,nnz(fnode2vox),0,0); % nothing
1076                 rnode2tet = test_vnode_in_tet(vox,tet);
1077                 unit_test_cmp(txt,nnz(rnode2tet),0,0); % nothing
1078                 mk_grid_c2f(tet,vox);
1079 
1080             case 5 % tet_edge on vox_face
1081                 txt = 'tet_edge on vox_face';
1082                 tet.nodes(:,1:2) = tet.nodes(:,1:2) - 0.3;
1083                 subplot(X,Y,i), show_test(vox,tet);
1084                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1085                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1086                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1087                 unit_test_cmp(txt,size(unique(intpts,'rows'),1),3,0); % 3 points
1088                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1089                 unit_test_cmp(txt,size(intpts,1),2,0); % 2 points
1090                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1091                 unit_test_cmp(txt,nnz(fnode2vox),0,0); % nothing
1092                 rnode2tet = test_vnode_in_tet(vox,tet);
1093                 unit_test_cmp(txt,nnz(rnode2tet),1,0); % 1 point
1094                 mk_grid_c2f(tet,vox);
1095 
1096             case 6 
1097                 txt = 'vox_node on tet_face';
1098                 tet.nodes(:,1:2) = tet.nodes(:,1:2) - 0.3;
1099                 tet.nodes(:,3)   = tet.nodes(:,3) -0.4;
1100                 subplot(X,Y,i), show_test(vox,tet);
1101                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1102                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1103                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1104                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1105                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1106                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1107                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1108                 unit_test_cmp(txt,nnz(fnode2vox),0,0); % nothing
1109                 rnode2tet = test_vnode_in_tet(vox,tet);
1110                 unit_test_cmp(txt,nnz(rnode2tet),1,0); % 1 point
1111                 mk_grid_c2f(tet,vox);
1112 
1113             case 7
1114                 txt = 'tet_node on vox_face';
1115                 tet.nodes(:,1) = tet.nodes(:,1) - 1;
1116                 tet.nodes(:,2:3)  = tet.nodes(:,2:3) + .5;
1117                 subplot(X,Y,i), show_test(vox,tet);
1118                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1119                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1120                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1121                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1122                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1123                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1124                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1125                 unit_test_cmp(txt,nnz(fnode2vox),1,0); % 1 point
1126                 rnode2tet = test_vnode_in_tet(vox,tet);
1127                 unit_test_cmp(txt,nnz(rnode2tet),0,0); % nothing
1128                 mk_grid_c2f(tet,vox);
1129 
1130 
1131             case 8
1132                 txt = 'tet_node on vox_edge';
1133                 tet.nodes(:,1) = tet.nodes(:,1) - 1;
1134                 tet.nodes(:,2)  = tet.nodes(:,2) + .5;
1135                 subplot(X,Y,i), show_test(vox,tet);
1136                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1137                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1138                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1139                 unit_test_cmp(txt,size(uniquetol(intpts,eps,'ByRows',true),1),1,0); % 1 point
1140                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1141                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1142                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1143                 unit_test_cmp(txt,nnz(fnode2vox),1,0); % 1 point
1144                 rnode2tet = test_vnode_in_tet(vox,tet);
1145                 unit_test_cmp(txt,nnz(rnode2tet),0,0); % nothing
1146                 mk_grid_c2f(tet,vox);
1147 
1148             case 9
1149                 txt = 'all nodes common';
1150                 subplot(X,Y,i), show_test(vox,tet);
1151                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1152                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1153                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1154                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1155                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1156                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1157                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1158                 unit_test_cmp(txt,nnz(fnode2vox),4,0); % 4 points
1159                 rnode2tet = test_vnode_in_tet(vox,tet);
1160                 unit_test_cmp(txt,nnz(rnode2tet),0,0); % nothing
1161                 mk_grid_c2f(tet,vox);
1162 
1163             case 10
1164                 txt = 'tet in vox';
1165                 tet.nodes = .25 + .5*tet.nodes;
1166                 subplot(X,Y,i), show_test(vox,tet);
1167                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1168                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1169                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1170                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1171                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1172                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1173                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1174                 unit_test_cmp(txt,nnz(fnode2vox),4,0); % 4 points
1175                 rnode2tet = test_vnode_in_tet(vox,tet);
1176                 unit_test_cmp(txt,nnz(rnode2tet),0,0); % nothing
1177                 mk_grid_c2f(tet,vox);
1178                 c2f = mk_grid_c2f(tet,vox);
1179                 max_vox = max(vox.nodes);
1180                 min_vox = min(vox.nodes);
1181                 edges = max_vox-min_vox;
1182                 vox_vol = edges(:,1) .* edges(:,2) .* edges(:,3);
1183                 tet_vol = get_elem_volume(tet);
1184                 unit_test_cmp(txt,c2f, 1, 0);
1185 
1186             case 11
1187                 txt = 'vox in tet';
1188                 tet.nodes =  4*tet.nodes - .25;
1189                 subplot(X,Y,i), show_test(vox,tet);
1190                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1191                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1192                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1193                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1194                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1195                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1196                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1197                 unit_test_cmp(txt,nnz(fnode2vox),0,0);  % nothing
1198                 rnode2tet = test_vnode_in_tet(vox,tet);
1199                 unit_test_cmp(txt,nnz(rnode2tet),8,0); % 8 points
1200                 c2f = mk_grid_c2f(tet,vox);
1201                 max_vox = max(vox.nodes);
1202                 min_vox = min(vox.nodes);
1203                 edges = max_vox-min_vox;
1204                 vox_vol = edges(:,1) .* edges(:,2) .* edges(:,3);
1205                 tet_vol = get_elem_volume(tet);
1206                 unit_test_cmp(txt,c2f, vox_vol / tet_vol, 0);
1207 
1208 
1209             case 12 
1210                 txt = 'tet_edge v. vox_face';
1211                 tet.nodes(:,1:2) = tet.nodes(:,1:2) - 0.3;
1212                 tet.nodes(:,3)   = tet.nodes(:,3) + .4;
1213                 subplot(X,Y,i), show_test(vox,tet);
1214                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1215                 unit_test_cmp(txt,size(intpts,1),2,0); % 2 points
1216                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1217                 unit_test_cmp(txt,size(intpts,1),2,0); % 2 points
1218                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1219                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1220                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1221                 unit_test_cmp(txt,nnz(fnode2vox),0,0); % nothing
1222                 rnode2tet = test_vnode_in_tet(vox,tet);
1223                 unit_test_cmp(txt,nnz(rnode2tet),0,0); % nothing
1224                 mk_grid_c2f(tet,vox);
1225 
1226 
1227             case 13
1228                 txt = 'everything';
1229                 tet.nodes = tet.nodes + 0.7;
1230                 subplot(X,Y,i), show_test(vox,tet);
1231                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1232                 unit_test_cmp(txt,size(intpts,1),3,0); % 3 points
1233                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1234                 unit_test_cmp(txt,size(intpts,1),3,0); % 3 points
1235                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1236                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1237                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1238                 unit_test_cmp(txt,nnz(fnode2vox),1,0); % 1 point
1239                 rnode2tet = test_vnode_in_tet(vox,tet);
1240                 unit_test_cmp(txt,nnz(rnode2tet),1,0); % 1 point
1241                 mk_grid_c2f(tet,vox);
1242 
1243                 %------- degenerate cases-----------%
1244             case 14 % common edge
1245                 txt = 'DG1: Common edge';
1246                 tet.nodes(:,1:2) = tet.nodes(:,1:2) + 1;
1247                 subplot(X,Y,i), show_test(vox,tet);
1248                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1249                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1250                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1251                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1252                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1253                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1254                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1255                 unit_test_cmp(txt,fnode2vox,[1;0;0;1],0); % 2 points
1256                 rnode2tet = test_vnode_in_tet(vox,tet);
1257                 unit_test_cmp(txt,nnz(rnode2tet),0,0);  % nothing
1258                 mk_grid_c2f(tet,vox);
1259 
1260             case 15 % common edge
1261                 txt = 'DG2: Common edge';
1262                 tet.nodes(:,1:2) = tet.nodes(:,1:2) + 1;
1263                 tet.nodes(:,3) = tet.nodes(:,3) + 0.5;
1264                 subplot(X,Y,i), show_test(vox,tet);
1265                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1266                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1267                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1268                 unit_test_cmp(txt,size(unique(intpts,'rows'),1),1,1); % 1
1269                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1270                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1271                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1272                 unit_test_cmp(txt,nnz(fnode2vox),1,0); % 1 point
1273                 rnode2tet = test_vnode_in_tet(vox,tet);
1274                 unit_test_cmp(txt,nnz(rnode2tet),1,0);  % 1 point
1275                 mk_grid_c2f(tet,vox);
1276 
1277              case 16 % common edge
1278                 txt = 'DG3: Common edge';
1279                 tet.nodes(:,1:2) = tet.nodes(:,1:2) + 1;
1280                 z = tet.nodes(:,3) == 0;
1281                 tet.nodes(z,3) = -.5;
1282                 tet.nodes(~z,3) = 1.5;
1283                 subplot(X,Y,i), show_test(vox,tet);
1284                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1285                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1286                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1287                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1288                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1289                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1290                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1291                 unit_test_cmp(txt,nnz(fnode2vox),0,0); % nothing
1292                 rnode2tet = test_vnode_in_tet(vox,tet);
1293                 unit_test_cmp(txt,nnz(rnode2tet),2,0);  % 2 points
1294                 mk_grid_c2f(tet,vox);
1295 
1296             case 17 % edge on face
1297                 txt = 'DG4: edge on face';
1298                 tet.nodes = [0 0 0; 1 .5 0; 0 1 0; 1 .5 1];
1299                 tet.nodes(:,1) = tet.nodes(:,1) - 1;
1300                 subplot(X,Y,i), show_test(vox,tet);
1301                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1302                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1303                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1304                 unit_test_cmp(txt,size(unique(intpts,'rows'),1),2,0); % 2
1305                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1306                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1307                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1308                 unit_test_cmp(txt,nnz(fnode2vox),2,0); % 2 points
1309                 rnode2tet = test_vnode_in_tet(vox,tet);
1310                 unit_test_cmp(txt,nnz(rnode2tet),0,0);  % nothing
1311                 mk_grid_c2f(tet,vox);
1312 
1313             case 18 % edge2edge only
1314                 txt = 'DG5: edge2edge only';
1315                 tet.nodes = [0 0 0; 1 .5 0; 0 1 0; 1 .5 1];
1316                 tet.nodes(:,1) = tet.nodes(:,1) - 1;
1317                 z = tet.nodes(:,3) == 0;
1318                 tet.nodes(z,3) = -.5;
1319                 tet.nodes(~z,3) = 1.5;
1320                 subplot(X,Y,i), show_test(vox,tet);
1321                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1322                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1323                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1324                 unit_test_cmp(txt,size(unique(intpts,'rows'),1),2,0); % 2 points
1325                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1326                 unit_test_cmp(txt,size(intpts,1),2,0); % 2 points
1327                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1328                 unit_test_cmp(txt,nnz(fnode2vox),0,0); % nothing
1329                 rnode2tet = test_vnode_in_tet(vox,tet);
1330                 unit_test_cmp(txt,nnz(rnode2tet),0,0);  % nothing
1331                 mk_grid_c2f(tet,vox);
1332 
1333             case 19 % edge on face
1334                 txt = 'DG6: edge on face';
1335                 tet.nodes = [0 0 0; 1 .5 0; 0 1 0; 1 .5 1];
1336                 tet.nodes(:,1) = tet.nodes(:,1) - 1;
1337                 tet.nodes(:,3) = tet.nodes(:,3) - .5;
1338                 subplot(X,Y,i), show_test(vox,tet);
1339                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1340                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1341                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1342                 unit_test_cmp(txt,size(unique(intpts,'rows'),1),1,0); % 1 point
1343                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1344                 unit_test_cmp(txt,size(intpts,1),1,0); % 1 point
1345                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1346                 unit_test_cmp(txt,nnz(fnode2vox),1,0); % 1 points
1347                 rnode2tet = test_vnode_in_tet(vox,tet);
1348                 unit_test_cmp(txt,nnz(rnode2tet),0,0);  % nothing
1349                 mk_grid_c2f(tet,vox);
1350 
1351             case 20 %face on face
1352                 txt = 'DG7: face on face';
1353                 tet.nodes(:,1) = tet.nodes(:,1) + 1;
1354                 tet.nodes(:,2) = tet.nodes(:,2) + 0.5;
1355                 tet.nodes(:,3) = tet.nodes(:,3) + 0.5;
1356                 subplot(X,Y,i), show_test(vox,tet);
1357                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1358                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1359                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1360                 unit_test_cmp(txt,size(uniquetol(intpts,2*eps,'ByRows',true),1),3,0); % 3 points
1361                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1362                 unit_test_cmp(txt,size(unique(intpts,'rows'),1),2,0); % 2 unique points (these edges are counted 3 times, but vox2edge takes care of this)
1363                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1364                 unit_test_cmp(txt,nnz(fnode2vox),1,0); % 1 point
1365                 rnode2tet = test_vnode_in_tet(vox,tet);
1366                 unit_test_cmp(txt,nnz(rnode2tet),1,0);  % 1 point
1367                 mk_grid_c2f(tet,vox);
1368 
1369            case 21 %face on face
1370                 txt = 'DG8: face on face';
1371                 tet.nodes(:,1) = tet.nodes(:,1) + 1;
1372                 tet.nodes(:,2) = tet.nodes(:,2) + 0.4;
1373                 tet.nodes(:,3) = tet.nodes(:,3) - 0.6;
1374                 subplot(X,Y,i), show_test(vox,tet);
1375                 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1376                 unit_test_cmp(txt,size(intpts,1),0,0); % nothing
1377                 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1378                 unit_test_cmp(txt,size(unique(intpts,'rows'),1),2,0); % 2 points
1379                 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1380                 unit_test_cmp(txt,size(unique(intpts,'rows'),1),2,0); % 2 unique points (these edges are counted 3 times, but vox2edge takes care of this)
1381                 [fnode2vox] = test_tnode_in_vox(vox,tet);
1382                 unit_test_cmp(txt,nnz(fnode2vox),1,0); % 1 point
1383                 rnode2tet = test_vnode_in_tet(vox,tet);
1384                 unit_test_cmp(txt,nnz(rnode2tet),0,0);  % nothing
1385                 mk_grid_c2f(tet,vox);
1386 
1387             otherwise
1388                 break;
1389         end
1390     end
1391     eidors_msg('log_level',ll);
1392 
1393 function do_small_test
1394    fmdl= ng_mk_cyl_models([2,2,.1],[16,1],[.1,0,.025]);
1395    xvec = [-1.5:1:1.5];
1396    yvec = [-1.6:.8:1.6];
1397    zvec = -1:1:2;
1398    rmdl = mk_grid_model([],xvec,yvec,zvec);
1399 
1400    eidors_cache clear
1401    tic
1402    c2f_a = mk_grid_c2f(fmdl, rmdl);
1403    toc
1404    eidors_cache clear
1405    tic
1406    opt.save_memory = 2;
1407    c2f_b = mk_grid_c2f(fmdl, rmdl, opt);
1408    toc
1409    eidors_cache clear
1410    tic
1411    c2f_c = mk_approx_c2f(fmdl, rmdl);
1412    toc
1413    assert(any(c2f_a(:) ~= c2f_b(:)) == 0);
1414    % the old function cannot be tested this way, it doesn't give the right
1415    % answer
1416    % assert(any(c2f_a(:) ~= c2f_c(:)) == 0);
1417    
1418 function do_realistic_test
1419 fmdl= ng_mk_cyl_models([2,2,.1],[16,1],[.1,0,.025]);
1420 xvec = [-1.5 -.5:.2:.5 1.5];
1421 yvec = [-1.6 -1:.2:1 1.6];
1422 zvec = 0:.25:2;
1423 rmdl = mk_grid_model([],xvec,yvec,zvec);
1424 
1425 tic
1426 opt.save_memory = 0;
1427 c2f_a = mk_grid_c2f(fmdl, rmdl,opt);
1428 t = toc;
1429 fprintf('Analytic: t=%f s\n',t);
1430 
1431 tic
1432 c2f_n = mk_approx_c2f(fmdl,rmdl);
1433 t = toc;
1434 fprintf('Approximate: t=%f s\n',t);
1435 
1436 
1437 tetvol = get_elem_volume(fmdl);
1438 opt = parse_opts(fmdl,rmdl);
1439 m = prepare_vox_mdl(rmdl,opt);
1440 
1441 
1442 f2c_a = bsxfun(@times, c2f_a, tetvol);
1443 f2c_a = bsxfun(@rdivide,f2c_a', m.volume); 
1444 img = mk_image(rmdl,0);
1445 img.elem_data = f2c_a*ones(size(fmdl.elems,1),1);
1446 subplot(132)
1447 show_fem(img);
1448 
1449 f2c_n = bsxfun(@times, c2f_n, tetvol);
1450 f2c_n = bsxfun(@rdivide,f2c_n', m.volume); 
1451 img = mk_image(rmdl,0);
1452 img.elem_data = f2c_n*ones(size(fmdl.elems,1),1);
1453 subplot(133)
1454 show_fem(img);
1455 subplot(131);
1456 h = show_fem(fmdl);
1457 set(h,'LineWidth',0.1)
1458 hold on
1459 h = show_fem(rmdl);
1460 set(h,'EdgeColor','b','LineWidth',2);
1461 hold off
1462 
1463 function do_edge2edge_timing_test
1464     fmdl= ng_mk_cyl_models([2,2],[8,1],[.15,0,.05]);
1465     xvec = linspace(-2,2,33);
1466     yvec = linspace(-2,2,33);
1467     zvec = 0:.5:2;
1468     rmdl = mk_grid_model([],xvec,yvec,zvec);
1469     tic
1470     test_tedge2vedge_intersections(rmdl,fmdl);
1471     toc
1472     
1473 function rnode2tet = test_vnode_in_tet(rmdl,fmdl)
1474     opt  = parse_opts(fmdl,rmdl);
1475     fmdl = prepare_fmdl(fmdl);
1476     rnode2tet = get_nodes_in_tets(fmdl,rmdl, opt);
1477 
1478 function [fnode2vox] = test_tnode_in_vox(rmdl,fmdl)
1479     fmdl = prepare_fmdl(fmdl);
1480     [fnode2vox] = get_nodes_in_voxels(fmdl,rmdl);
1481     
1482 function [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(rmdl,fmdl)
1483     opt = parse_opts(fmdl,rmdl);
1484     m = prepare_vox_mdl(rmdl,opt);
1485     fmdl = prepare_fmdl(fmdl);
1486     [intpts, face2edge, face2intpt, edge2intpt] = get_voxel_intersection_points(fmdl,m.faces,opt);
1487         
1488 function [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(rmdl,fmdl)
1489     opt = parse_opts(fmdl,rmdl);
1490     m = prepare_vox_mdl(rmdl,opt);
1491     fmdl = prepare_fmdl(fmdl);
1492     [intpts, face2edge, face2intpt, edge2intpt] = get_tet_intersection_points(fmdl,m,opt);
1493     
1494 function [intpts, tedge2vedge, tedge2intpt, vedge2intpt] = test_tedge2vedge_intersections(rmdl,fmdl)
1495     opt = parse_opts(fmdl,rmdl);
1496     m = prepare_vox_mdl(rmdl,opt);
1497     fmdl = prepare_fmdl(fmdl);
1498     [intpts, tedge2vedge, tedge2intpt, vedge2intpt] = find_edge2edge_intersections(fmdl.edges,fmdl.nodes,m.edges,rmdl.nodes, opt.tol_node2tet);
1499 
1500 function show_test(vox,tet)
1501     show_fem(vox);
1502     hold on
1503     h = show_fem(tet);
1504     set(h, 'EdgeColor','b');
1505     hold off
1506     axis auto