function [t]=MyCrust(p)
%% Main
starttime=clock;
%add points to the given ones, this is usefull to create outside tetraedrom
tic
p=AddShield(p);
fprintf('Addedded Shield: %4.4f s\n',toc)
tic
tetr=delaunayn(p);%creating tedraedron
tetr=int32(tetr);%use integer to save memory
fprintf('Delaunay Triangulation Time: %4.4f s\n',toc)
%connectivity data
%find triangles to tetraedrom and tetraedrom to triangles connectivity data
tic
[t2tetr,tetr2t]=Connectivity(tetr);
fprintf('Connectivity Time: %4.4f s\n',toc)
tic
[cc,r]=CC();%Circumcenters of tetraedroms
fprintf('Circumcenters Time: %4.4f s\n',toc)
clear n
tic
t=Walking();%Flagging tetraedroms as inside or outside
fprintf('Walking Time: %4.4f s\n',toc)
time=etime(clock,starttime);
fprintf('Total Time: %4.4f s\n',time)
%% Circumcenters(Nested嵌套)
function [cc,r]=CC()
%finds circumcenters fro a set of tetraedrom
%points of tetraedrom
p1=(p(tetr(:,1),:));
p2=(p(tetr(:,2),:));
p3=(p(tetr(:,3),:));
p4=(p(tetr(:,4),:));
%vectors of tetraedrom edges
v21=p(tetr(:,1),:)-p(tetr(:,2),:);
v31=p(tetr(:,3),:)-p(tetr(:,1),:);
v41=p(tetr(:,4),:)-p(tetr(:,1),:);
%preallocation
cc=zeros(size(tetr,1),3);
%Solve the system using cramer method
d1=sum(v41.*(p1+p4)*.5,2);
d2=sum(v21.*(p1+p2)*.5,2);
d3=sum(v31.*(p1+p3)*.5,2);
det23=(v21(:,2).*v31(:,3))-(v21(:,3).*v31(:,2));
det13=(v21(:,3).*v31(:,1))-(v21(:,1).*v31(:,3));
det12=(v21(:,1).*v31(:,2))-(v21(:,2).*v31(:,1));
Det=v41(:,1).*det23+v41(:,2).*det13+v41(:,3).*det12;
detx=d1.*det23+...
v41(:,2).*(-(d2.*v31(:,3))+(v21(:,3).*d3))+...
v41(:,3).*((d2.*v31(:,2))-(v21(:,2).*d3));
dety=v41(:,1).*((d2.*v31(:,3))-(v21(:,3).*d3))+...
d1.*det13+...
v41(:,3).*((d3.*v21(:,1))-(v31(:,1).*d2));
detz=v41(:,1).*((v21(:,2).*d3)-(d2.*v31(:,2)))...
+v41(:,2).*(d2.*v31(:,1)-v21(:,1).*d3)...
+d1.*(det12);
%Circumcenters
cc(:,1)=detx./Det;
cc(:,2)=dety./Det;
cc(:,3)=detz./Det;
%Circumradius
r=((sum((p2-cc).^2,2))).^.5;%quadrato del raggio
end
%% Connectivity(Nested)
function [t2tetr,tetr2t]=Connectivity(tetr)
numt = size(tetr,1);
vect = 1:numt;
t = [tetr(:,[1,2,3]); tetr(:,[2,3,4]); tetr(:,[1,3,4]);tetr(:,[1,2,4])];%triangles not unique
[t,j,j] = unique(sort(t,2),'rows');%triangles
t2tetr = [j(vect), j(vect+numt), j(vect+2*numt),j(vect+3*numt)];%each tetraedrom has 3 triangles
% triang-to-tetr connectivity
nume = size(t,1);
tetr2t = zeros(nume,2,'int32');
count= ones(nume,1,'int8');
for k = 1:numt
for j=1:4
ce = t2tetr(k,j);
tetr2t(ce,count(ce)) = k;
count(ce)=count(ce)+1;
end
end
end % connectivity()
%% Walking(Nested)
function t=Walking()
np=size(p,1)-540;%540 = number of shield points put at the end of array 放在数组后面的保护点数目
numtetr=size(tetr,1);
nt=size(tetr2t,1);
deleted=true(numtetr,1);%deleted tetraedroms
checked=false(numtetr,1);%checked tetraedros
onfront=false(nt,1);%tetraedroms that need to be checked
countchecked=0;
%First flag as outsde tetraedroms with Shield points
for i=1:numtetr
for j=1:4
if tetr(i,j)>np;
deleted(i)=true;
checked(i)=true;
onfront(t2tetr(i,:))=true;
countchecked=countchecked+1;
break
end
end
end
%tollerances to mark as in or out
toll=zeros(nt,1)+.95;
level=0;
alpha=zeros(nt,1);%intersection factor
%it is computed from radius of the tetraedroms circumscribed sphere
% and the distance between their center
for i=1:nt
if tetr2t(i,2)>0 %jump boundary tetraedrom
distcc=sum((cc(tetr2t(i,1),:)-cc(tetr2t(i,2),:)).^2,2);%distance from circumcenters
%intersection factor
alpha(i)=(-distcc+r(tetr2t(i,1))^2+r(tetr2t(i,2))^2)/(2*r(tetr2t(i,1))*r(tetr2t(i,2)));
end
end
clear cc
% Now we scan扫描 all tetraedroms. When one is scanned put on front is
% neighobur. This means that now even the neighobour can be
% checked. At the begining only tetraedroms with shield points are
% on front, because we are sure the are out.
% Tetraedrom with high intersction factor will be marked as equal
% else different.
% When I say high i mean under a set tollerance that becames lower as the algorithm
% progresses. This Aims to avoid errors propagation when a tetraedrom is wrong marked.
%
while countchecked<numtetr && level<50
level=level+1;%level of scan reached
for id=1:nt%loop trough triangles
if onfront(id)
tetr1=tetr2t(id,1);tetr2=tetr2t(id,2);%tetraedroms linked to triangle under analysis
if tetr2==0 %do not check boundary triangles
onfront(id)=false;
continue
elseif (checked(tetr1) && checked(tetr2)) %tetraedroms are already checked
onfront(id)=false;
continue
end
if alpha(id)>=toll(id) %flag as equal
if checked(tetr1)%find the checked one between the two
deleted(tetr2)=deleted(tetr1) ;%flag as equal
checked(tetr2)=true;%check
countchecked=countchecked+1;
onfront(t2tetr(tetr2,:))=true;%put on front all tetreadrom triangles
else
deleted(tetr1)=deleted(tetr2) ;%flag as equal
checked(tetr1)=true;%check
countchecked=countchecked+1;
onfront(t2tetr(tetr1,:))=true;%put on front all tetreadrom triangles
end
onfront(id)=false;%remove from front
elseif alpha(id)<-toll(id)%flag as different
if checked(tetr1)%find the checked one between the two
deleted(tetr2)=~(deleted(tetr1)) ;%flag as different
checked(tetr2)=true;%check
countchecked=countchecked+1;
onfront(t2tetr(tetr2,:))=true;%put on front all tetreadrom triangles
else
deleted(tetr1)=~(deleted(tetr2)) ;%flag as different
checked(tetr1)=true;%check
countchecked=countchecked+1;
onfront(t2tetr(tetr1,:))=true;%put on front all tetreadrom triangles
end
onfront(id)=false;%remove from front
else
toll(id)=toll(id)-.05;%tolleraces were too high next time will be lower
end
end
end
if l
