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 6423 2022-11-30 11:44:48Z 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 % Octave <=7.3 gives massive text output
0043 % from Qhull. Instead go straight to joggle
0044 % version
0045 
0046 % TODO: should we always use joggle anyway?
0047 if exist('OCTAVE_VERSION')
0048     [K, V] = convhulln(pts,{'Qt Pp Qs QJ'});
0049     if V<1e-8; V=0; end
0050 else
0051     try
0052         [K, V] = convhulln(pts,{'Qt Pp Qs'});
0053     catch
0054         %redo it with "Joggle", but set to zero if small
0055         [K, V] = convhulln(pts,{'Qt Pp Qs QJ'});
0056         if V<1e-8; V=0; end
0057     end
0058 end
0059    
0060 
0061 function [K,V] = call_convhulln_old(pts,p);
0062   K=0; V=0;
0063   dim = size(pts,2);
0064 
0065 try
0066        [K, V] = convhulln(pts,{'Qt Pp Qs'});
0067 catch err
0068   ok = false;
0069   if exist('OCTAVE_VERSION')
0070      if strcmp(err.message,'convhulln: qhull failed')
0071         err.identifier =  'MATLAB:qhullmx:DegenerateData';
0072      end
0073         
0074   end
0075   switch err.identifier
0076      case {'MATLAB:qhullmx:DegenerateData', 'MATLAB:qhullmx:UndefinedError'}
0077         % border case point may be included multiple times.
0078         % this is OK... though we may miss cases where more
0079         % points should have been found but were not
0080         u = uniquetol(pts*p.scale,1e-14,'ByRows',true,'DataScale', 1);
0081         ok = ok | (size(u,1) <= dim  );
0082         if ~ok; switch dim;
0083            case 2; ok = colinear_test2d(u);
0084            case 3; ok = colinear_test3d(pts*p.scale);
0085            otherwise; error('not 2D or 3D');
0086         end; end
0087   end
0088 %    Save cases were errors called
0089 %      load -mat CHP.mat ptsi;
0090 %      ptsi{end+1} = pts;
0091 %      save -mat CHP.mat ptsi;
0092   if ~ok
0093      if      eidors_debug('query','mk_tet_c2f:convhulln');
0094         debug_plot_tet(p.fmdl,p.rmdl,p.tri_todo,p.t, p.pts)
0095         keyboard
0096      elseif  eidors_debug('query','mk_tri2tet_c2f:convhulln');
0097         debug_plot_tri2tet(p.fmdl,p.rmdl,p.v,p.t, p.bot, p.top, p.pts)
0098         keyboard
0099      else
0100         fprintf('\n');
0101         eidors_msg(['convhulln has thrown an error. (',err.message,')', ...
0102            'Enable "eidors_debug on convhulln_clean" and re-run to see a debug plot'],0);
0103         rethrow(err);
0104      end
0105   end
0106 end
0107 
0108 % test for colinearity in 2D
0109 function ok = colinear_test2d(u,ok) 
0110    ok = false;
0111    cp = bsxfun(@minus, u(2:end,:), u(1,:));
0112    l = sqrt(sum(cp.^2,2));
0113    cp = bsxfun(@rdivide, cp, l);
0114    u = uniquetol(cp,1e-14,'ByRows',true,'DataScale',1);
0115    ok = ok | size(u,1) == 1; % co-linear points
0116 
0117 % test for colinearity in 3D
0118 function ok = colinear_test3d(pts);
0119    ok = false;
0120    u12 = uniquetol(pts(:,1:2),1e-14,'ByRows',true,'DataScale',1);
0121    cp = bsxfun(@minus, u12(2:end,1:2), u12(1,1:2));
0122    l = sqrt(sum(cp.^2,2));
0123    cp = bsxfun(@rdivide, cp, l);
0124    % counteract colinear vectors in different directions
0125    cp = abs(cp); 
0126    un = uniquetol(cp,1e-12,'ByRows',true,'DataScale',1);
0127    ok = ok | size(un,1) == 1; % co-linear points
0128    if ok; return; end
0129 
0130    % test if all points lie on the top or bottom caps
0131    top = max(pts(:,3));
0132    bot = min(pts(:,3));
0133    ok = ok | all(abs(pts(:,3) - top) < eps);
0134    ok = ok | all(abs(pts(:,3) - bot) < eps);
0135 
0136 function debug_plot_tet(fmdl,rmdl,tri_todo,t, pts)
0137    clf
0138    tri.nodes = fmdl.nodes;
0139    vox.nodes = rmdl.nodes;
0140    tri.type = 'fwd_model';
0141    vox.type = 'fwd_model';
0142    vox.elems = rmdl.elems(v,:);
0143    vox.boundary = p.vox.elems;
0144    tri.elems = fmdl.elems(tri_todo(t),:);
0145    show_fem(vox)
0146    hold on
0147    h = show_fem(tri);
0148    set(h,'EdgeColor','b')
0149    pts = bsxfun(@plus,pts*scale,ctr);
0150    plot(pts(:,1),pts(:,2),'o');
0151    hold off
0152    axis auto
0153 
0154 function debug_plot_tri2tet(fmdl,rmdl,v,t, bot, top, pts)
0155    clf
0156    tet.nodes = fmdl.nodes;
0157    tri.nodes = repmat(rmdl.nodes(rmdl.elems(v,:),:),2,1);
0158    tri.nodes(:,3) = [repmat(bot,3,1); repmat(top,3,1)];
0159    tri.elems = [ 1 2 5 4
0160                  2 3 6 5
0161                  3 1 4 6];
0162    tri.boundary = tri.elems;
0163    tet.type = 'fwd_model';
0164    tri.type = 'fwd_model';
0165    tet.elems = fmdl.elems(t,:);
0166    clf
0167    show_fem(tri)
0168    hold on
0169    h = show_fem(tet);
0170    set(h,'EdgeColor','b')
0171 %    pts = bsxfun(@plus,pts*scale,ctr);
0172    plot3(pts(:,1),pts(:,2),pts(:,3),'o');
0173    hold off
0174    axis auto
0175 
0176 function do_unit_test
0177    t1.type = 'fwd_model'; t1.elems = [1 2 3 4];
0178    t1.nodes = [0 0 0; 0 1 0; 1 0 0; 0 0 1]; t2 = t1;
0179    unit_test_cmp('A',mk_tet_c2f(t1,t2),    1,10*eps);
0180    t2.nodes(end,end) = -1;
0181    unit_test_cmp('B',mk_tet_c2f(t1,t2),    0,10*eps);
0182