convhulln_clean

 CONVHULLN_CLEAN: run convhulln and catch errors

0001 function [K,V] = convhulln_clean(pts,p);
0002 % CONVHULLN_CLEAN: run convhulln and catch errors
0003 
0004 % (C) 2018 Bartlomiej Grychtol, Andy Adler
0005 % License: GPL version 2 or 3
0006 % $Id: convhulln_clean.m 5890 2018-12-23 21:55:26Z aadler $
0007 
0008   if ischar(pts) && strcmp(pts,'UNIT_TEST'); do_unit_test; return; end
0009 
0010 
0011   K=0; V=0; dim = size(pts,2);
0012 
0013   if size(pts,1) < 3; return; end 
0014   % move points to origin (helps for small elements at
0015   % large coordinates
0016   ctr = mean(pts);
0017   pts = bsxfun(@minus,pts,ctr);
0018   scale = max(abs(pts(:)));
0019 
0020   if scale == 0; return; end  %when there's only one point
0021 
0022   % scale largest coordinate to 1 (helps with precision)
0023   pts = pts ./ scale;
0024   p.scale = scale;
0025 
0026   % force thorough search for initinal simplex and
0027   % supress precision warnings
0028   pts= uniquetol(pts,1e-13,'ByRows',true,'DataScale',1);
0029   % This won't catch all cases
0030   if size(pts,1)<=size(pts,2); return; end
0031   if any(std(pts)<1e-14); return; end
0032   [K,V] = call_convhulln(pts,p);
0033 
0034  % numerical issues may produce tiny negative volume
0035  % undo scaling
0036   V = scale^dim * max(V,0);
0037 
0038 function [K,V] = call_convhulln(pts,p);
0039   K=0; V=0;
0040   dim = size(pts,2);
0041 
0042 try
0043     [K, V] = convhulln(pts,{'Qt Pp Qs'});
0044 catch
0045     %redo it with "Joggle", but set to zero if small
0046     [K, V] = convhulln(pts,{'Qt Pp Qs QJ'});
0047     if V<1e-8; V=0; end
0048 end
0049    
0050 
0051 function [K,V] = call_convhulln_old(pts,p);
0052   K=0; V=0;
0053   dim = size(pts,2);
0054 
0055 try
0056        [K, V] = convhulln(pts,{'Qt Pp Qs'});
0057 catch err
0058   ok = false;
0059   if exist('OCTAVE_VERSION')
0060      if strcmp(err.message,'convhulln: qhull failed')
0061         err.identifier =  'MATLAB:qhullmx:DegenerateData';
0062      end
0063         
0064   end
0065   switch err.identifier
0066      case {'MATLAB:qhullmx:DegenerateData', 'MATLAB:qhullmx:UndefinedError'}
0067         % border case point may be included multiple times.
0068         % this is OK... though we may miss cases where more
0069         % points should have been found but were not
0070         u = uniquetol(pts*p.scale,1e-14,'ByRows',true,'DataScale', 1);
0071         ok = ok | (size(u,1) <= dim  );
0072         if ~ok; switch dim;
0073            case 2; ok = colinear_test2d(u);
0074            case 3; ok = colinear_test3d(pts*p.scale);
0075            otherwise; error('not 2D or 3D');
0076         end; end
0077   end
0078 %    Save cases were errors called
0079 %      load -mat CHP.mat ptsi;
0080 %      ptsi{end+1} = pts;
0081 %      save -mat CHP.mat ptsi;
0082   if ~ok
0083      if      eidors_debug('query','mk_tet_c2f:convhulln');
0084         debug_plot_tet(p.fmdl,p.rmdl,p.tri_todo,p.t, p.pts)
0085         keyboard
0086      elseif  eidors_debug('query','mk_tri2tet_c2f:convhulln');
0087         debug_plot_tri2tet(p.fmdl,p.rmdl,p.v,p.t, p.bot, p.top, p.pts)
0088         keyboard
0089      else
0090         fprintf('\n');
0091         eidors_msg(['convhulln has thrown an error. (',err.message,')', ...
0092            'Enable "eidors_debug on convhulln_clean" and re-run to see a debug plot'],0);
0093         rethrow(err);
0094      end
0095   end
0096 end
0097 
0098 % test for colinearity in 2D
0099 function ok = colinear_test2d(u,ok) 
0100    ok = false;
0101    cp = bsxfun(@minus, u(2:end,:), u(1,:));
0102    l = sqrt(sum(cp.^2,2));
0103    cp = bsxfun(@rdivide, cp, l);
0104    u = uniquetol(cp,1e-14,'ByRows',true,'DataScale',1);
0105    ok = ok | size(u,1) == 1; % co-linear points
0106 
0107 % test for colinearity in 3D
0108 function ok = colinear_test3d(pts);
0109    ok = false;
0110    u12 = uniquetol(pts(:,1:2),1e-14,'ByRows',true,'DataScale',1);
0111    cp = bsxfun(@minus, u12(2:end,1:2), u12(1,1:2));
0112    l = sqrt(sum(cp.^2,2));
0113    cp = bsxfun(@rdivide, cp, l);
0114    % counteract colinear vectors in different directions
0115    cp = abs(cp); 
0116    un = uniquetol(cp,1e-12,'ByRows',true,'DataScale',1);
0117    ok = ok | size(un,1) == 1; % co-linear points
0118    if ok; return; end
0119 
0120    % test if all points lie on the top or bottom caps
0121    top = max(pts(:,3));
0122    bot = min(pts(:,3));
0123    ok = ok | all(abs(pts(:,3) - top) < eps);
0124    ok = ok | all(abs(pts(:,3) - bot) < eps);
0125 
0126 function debug_plot_tet(fmdl,rmdl,tri_todo,t, pts)
0127    clf
0128    tri.nodes = fmdl.nodes;
0129    vox.nodes = rmdl.nodes;
0130    tri.type = 'fwd_model';
0131    vox.type = 'fwd_model';
0132    vox.elems = rmdl.elems(v,:);
0133    vox.boundary = p.vox.elems;
0134    tri.elems = fmdl.elems(tri_todo(t),:);
0135    show_fem(vox)
0136    hold on
0137    h = show_fem(tri);
0138    set(h,'EdgeColor','b')
0139    pts = bsxfun(@plus,pts*scale,ctr);
0140    plot(pts(:,1),pts(:,2),'o');
0141    hold off
0142    axis auto
0143 
0144 function debug_plot_tri2tet(fmdl,rmdl,v,t, bot, top, pts)
0145    clf
0146    tet.nodes = fmdl.nodes;
0147    tri.nodes = repmat(rmdl.nodes(rmdl.elems(v,:),:),2,1);
0148    tri.nodes(:,3) = [repmat(bot,3,1); repmat(top,3,1)];
0149    tri.elems = [ 1 2 5 4
0150                  2 3 6 5
0151                  3 1 4 6];
0152    tri.boundary = tri.elems;
0153    tet.type = 'fwd_model';
0154    tri.type = 'fwd_model';
0155    tet.elems = fmdl.elems(t,:);
0156    clf
0157    show_fem(tri)
0158    hold on
0159    h = show_fem(tet);
0160    set(h,'EdgeColor','b')
0161 %    pts = bsxfun(@plus,pts*scale,ctr);
0162    plot3(pts(:,1),pts(:,2),pts(:,3),'o');
0163    hold off
0164    axis auto
0165 
0166 function do_unit_test
0167    t1.type = 'fwd_model'; t1.elems = [1 2 3 4];
0168    t1.nodes = [0 0 0; 0 1 0; 1 0 0; 0 0 1]; t2 = t1;
0169    unit_test_cmp('A',mk_tet_c2f(t1,t2),    1,10*eps);
0170    t2.nodes(end,end) = -1;
0171    unit_test_cmp('B',mk_tet_c2f(t1,t2),    0,10*eps);
0172