tet_to_inequal

 [A,b]=tet_to_inequal(v)
 Given the vertices of a simplex v return a system
 of linear inequalities so that a point x in in 
 the simplex iff Ax <= b

 [A,b]=tet_to_inequal(v)
 Given the Nx3 list of vertices v and the Nx4 list of elements e indexing
 v, return a system of linear inequalities so that a point x in in 
 the simplex iff Ax <= b

0001 function [A,b]=tet_to_inequal(v,e)
0002 % [A,b]=tet_to_inequal(v)
0003 % Given the vertices of a simplex v return a system
0004 % of linear inequalities so that a point x in in
0005 % the simplex iff Ax <= b
0006 %
0007 % [A,b]=tet_to_inequal(v)
0008 % Given the Nx3 list of vertices v and the Nx4 list of elements e indexing
0009 % v, return a system of linear inequalities so that a point x in in
0010 % the simplex iff Ax <= b
0011 
0012 % (C) 2012 Bill Lionheart. License GPL v2 or v3
0013 % $Id: tet_to_inequal.m 5112 2015-06-14 13:00:41Z aadler $
0014 
0015 if ischar(v) && strcmp(v,'UNIT_TEST'); do_unit_test; return; end
0016 
0017 if nargin == 1
0018    e = [1 2 3 4];
0019    ne = 1;
0020 end
0021 ne = size(e,1);
0022 d1 = size(e,2);
0023 d  = size(v,2);
0024 % [d1,d]=size(v(e,:));
0025 if d~=3
0026     error('Only works for dimension 3 at the moment')
0027 end
0028 if d1 ~= d+1
0029     error('Simplex needs to have one more vertex than dimension')
0030 end
0031 % the d edges at vertex 1
0032 edges1= v(e(:,2:end)',:) - kron(v(e(:,1),:), ones(3,1));%ones(3,1)*v(e(:,1),:);
0033 
0034 dt = my_det(edges1);
0035 idx = dt<0;
0036 if any(idx)
0037    e(idx, :) = [e(idx,2) e(idx,1) e(idx,3:end)];
0038    idx3 = reshape(repmat(idx,1,3)',1,[])';
0039    edges1(idx3,:)= v(e(idx,2:end)',:)-reshape(repmat(v(e(idx,1),:),1,3)',3,[])';
0040 end
0041 idx1 = 1:3*ne; % 1 2 3 4 5 6 ...
0042 idx3 = reshape(circshift(reshape(idx1,3,[]),-1),1,[]); % 2 3 1 5 6 4 ...
0043 idx2 = 1:4*ne; idx2(4:4:end) = [];
0044 A(idx2,:) = -reshape(my_cross(edges1(idx1,:),edges1(idx3,:)),[],3);
0045 b(idx2) = sum(A(idx2,:).*v(e(:,2:end)',:),2);
0046 % One more face not containing vertex 1
0047 edges3= [v(e(:,3)',:)-v(e(:,2)',:);v(e(:,4)',:)-v(e(:,2)',:)];
0048 idx4 = 4:4:4*ne;
0049 A(idx4,:) = reshape(my_cross(edges3(1:end/2,:),edges3(end/2+1:end,:)),[],3);
0050 b(idx4) = sum(A(idx4,:).*v(e(:,2)',:),2);
0051 b=b';
0052 end
0053 
0054 function c = my_cross(a,b)
0055 c = [a(:,2).*b(:,3)-a(:,3).*b(:,2) 
0056      a(:,3).*b(:,1)-a(:,1).*b(:,3)
0057      a(:,1).*b(:,2)-a(:,2).*b(:,1)];
0058 end  
0059 
0060 function d = my_det(a)
0061 ln = size(a,1);
0062 c1 = 1:3:ln;
0063 c2 = 2:3:ln;
0064 c3 = 3:3:ln;
0065 d = a(c1,1).*a(c2,2).*a(c3,3) + ...
0066     a(c2,1).*a(c3,2).*a(c1,3) + ...
0067     a(c3,1).*a(c1,2).*a(c2,3) - ...
0068     a(c3,1).*a(c2,2).*a(c1,3) - ...
0069     a(c2,1).*a(c1,2).*a(c3,3) - ...
0070     a(c1,1).*a(c3,2).*a(c2,3);
0071 end
0072 
0073 function do_unit_test
0074 simple_test
0075 test_3d
0076 end
0077 
0078 function test_3d
0079 imdl = mk_common_model('a3cr',16);
0080 fmdl = imdl.fwd_model;
0081 fmdl = linear_reorder(fmdl,1);
0082 [A1,b1] = tet_to_inequal(fmdl.nodes,fmdl.elems);
0083 v = fmdl.nodes(fmdl.elems(1,:),:);
0084 [A2,b2] = tet_to_inequal(v);
0085 unit_test_cmp('Parallel vs serial', A1((1:4),:),A2)
0086 unit_test_cmp('Parallel vs serial', b1(1:4),b2)
0087 v = fmdl.nodes(fmdl.elems(2,:),:);
0088 [A2,b2] = tet_to_inequal(v);
0089 unit_test_cmp('Parallel vs serial', A1(4*(2-1)+(1:4),:),A2)
0090 unit_test_cmp('Parallel vs serial', b1(4*(2-1)+(1:4)),b2)
0091 v = fmdl.nodes(fmdl.elems(57,:),:);
0092 [A2,b2] = tet_to_inequal(v);
0093 unit_test_cmp('Parallel vs serial', A1(4*(57-1)+(1:4),:),A2)
0094 unit_test_cmp('Parallel vs serial', b1(4*(57-1)+(1:4)),b2)
0095 end
0096 
0097 function simple_test
0098 v1=[0,0,0;eye(3)];
0099 [A1,b1] =  tet_to_inequal(v1);
0100 out = all( A1*[0.1;0.1;0.1]-b1<0);
0101 correct = 1;
0102 unit_test_cmp('Point in unit tetrahedron', out, correct)
0103 end