%%?MyCrust
%
%Simple?surface?recostruction?program?based?on?Crust?algorithm
%Given?a?set?of?3D?points?returns?a?triangulated?tight?surface.
%
%The?more?points?there?are?the?best?the?surface?will?be?fitted
%although?you?will?have?to?wait?more.?For?very?large?models?an?
%help?memory?errors?may?occurs.
%It?is?important?even?the?point?distribution?generally?uniformly?
%?distributed?points?with?denser?zones?in?high?curvature?features
%?give?the?best?results.
%?
%?Remember??crust?algorithom?needs?a?cloud?representing?a?volume?
%?so?open?surface?may?give?inaccurate?results.?
%??
%??
%?If?any?problems?occurs?in?execution?or?if?you?found?a?bug
%?have?a?suggestion?or?question?just?contact?me?at:??
%???
%?giaccariluigi@msn.com
%?
%
%
%
%Here?is?a?simple?example:
%
%load?Dino.mat%load?input?points?from?mat?file
%
%[t]=MyCrust（p）;
%
%?figure（1）
%?????????hold?on
%?????????title（‘Output?Triangulation‘‘fontsize‘14）
%?????????axis?equal
%?????????trisurf（tp（:1）p（:2）p（:3）‘facecolor‘‘c‘‘edgecolor‘‘b‘）
%
%Input:
%??????????????p?is?a?Nx3?array?containing?the?3D?set?of?points
%Output:
%??????????????t?are?points?id?contained?in?triangles?Nx3?array?too.
%
%?See?also?qhull?voronoin?convhulln?delaunay?delaunay3?tetramesh.
%
%Author:Giaccari?Luigi
%Last?Update:?1/12/2008
%Creation:?10/10/2008


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
[t2tetrtetr2t]=Connectivity（tetr）;
fprintf（‘Connectivity?Time:?%4.4f?s\n‘toc）
tic
[ccr]=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（clockstarttime）;
fprintf（‘Total?Time:?%4.4f?s\n‘time）





%%?Circumcenters（Nested）
????function?[ccr]=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（tetr1）3）;






?????????%Solve?the?system?using?cramer?method
????????d1=sum（v41.*（p1+p4）*.52）;
????????d2=sum（v21.*（p1+p2）*.52）;
????????d3=sum（v31.*（p1+p3）*.52）;

????????det23=（v21（:2）.*v31（:3））-（v21（:3）.*v31（:2））;
????????det13=（v21（:3）.*v31（:1））-（v21（:1）.*v31（:3））;
????????det12=（v2