mk_grid_c2f

PURPOSE ^

MK_GRID_C2F - calculate a coarse2fine mapping for grid coarse models.

SYNOPSIS ^

function [c2f, m] = mk_grid_c2f(fmdl, rmdl, opt)

DESCRIPTION ^

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

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