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