MATLAB 三维点云重构Surface recostruction from scattered points cloud

%% 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.  Surface with small holes are not
% considered open surface and generally are processed well. 
%
%
% 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 .
%
% See also qhull voronoin convhulln delaunay delaunay3 tetramesh.
%
%Author:Giaccari Luigi
%Last Update: 03/1/2008
%Created: 10/10/2008


function [t]=MyCrust（p）
%error check

if nargin>1
    error（‘The only input must be the Nx3 array of points‘）;
end

[n]=size（p2）;
if n ~=3
       error（‘Input 3D points must be stored in a Nx3 array‘）;
end 
clear  n


%%   Main
starttime=clock;

%add points to the given ones this is usefull
%to create outside tetraedroms
tic
[pnshield]=AddShield（p）;
fprintf（‘Added Shield: %4.4f s
‘toc）



tic
tetr=delaunayn（p）;%creating tedraedron
tetr=int32（tetr）;%use integer to save memory
fprintf（‘Delaunay Triangulation Time: %4.4f s
‘toc）



%connectivity data
%find triangles to tetraedrom and tetraedrom to triangles connectivity data
tic
[t2tetrtetr2tt]=Connectivity（tetr）;
fprintf（‘Connectivity Time: %4.4f s
‘toc）


tic
[ccr]=CC（ptetr）;%Circumcenters of tetraedroms
fprintf（‘Circumcenters Time: %4.4f s
‘toc）
clear n


tic
tbound=Marking（ptetrtetr2tt2tetrccrnshield）;%Flagging tetraedroms as inside or outside
fprintf（‘Walking Time: %4.4f s
‘toc）

%recostructed raw surface
t=t（tbound:）;

time=etime（clockstarttime）;
fprintf（‘Total Time: %4.4f s
‘time）


end


%% Circumcenters
function [ccr]=CC（ptetr）
%finds circumcenters from a set of tetraedroms

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

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2009-02-09 23:08  MyCrust090209
     文件       67998  2007-04-13 18:57  MyCrust090209Beethoven.mat
     文件       45713  2008-11-20 23:01  MyCrust090209Block.mat
     文件     1100582  2008-11-20 23:13  MyCrust090209Chair.mat
     文件      550432  2008-11-20 23:05  MyCrust090209Elephant.mat
     文件      290952  2008-07-02 00:24  MyCrust090209gargo50k.mat
     文件      188902  2008-06-28 20:52  MyCrust090209hippo.mat
     文件      985815  2008-11-20 23:13  MyCrust090209Horse.mat
     文件      218573  2009-01-05 18:37  MyCrust090209Knot.mat
     文件       12438  2009-02-09 22:48  MyCrust090209MyCrust.m
     文件      358743  2009-01-05 18:46  MyCrust090209OilPump.mat
     文件      306544  2008-11-20 23:11  MyCrust090209Skull.mat
     文件      804127  2008-11-20 23:10  MyCrust090209Standford_Bunny.mat
     文件         953  2009-02-09 23:14  MyCrust090209TestMyCrust.m
     文件     1121230  2009-01-22 23:46  MyCrust090209Vertebra1.mat
     文件        1335  2009-05-18 18:13  license.txt