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