Matlab 三维点云 三角化

%% 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