%% Reconstructing the 3D surface of  point cloud data
% The boundary or alphaShape function is used to construct the surface of the point cloud data, the subdivide function is used to refine the surface, and smoothSurfaceMesh is used to smooth the 3D surface.
% The output ply file can be imported from COMSOL geometry

clear
% load SEG_Saltxyz.mat
% trimesh(t,p(:,1),p(:,2),p(:,3))
% 
% ptCloud = pointCloud(p);

ptCloud = pcread('cave02m.ply');
figure ; pcshow(ptCloud,"MarkerSize",10)
 
%Cloud point downsampling
gridStep = 0.5;
ptCloudA = pcdownsample(ptCloud,'gridAverage',gridStep);
p = ptCloudA.Location;
p2 = double(p-mean(p));
% figure
% pcshow(ptCloudA,"MarkerSize",10);


%% Generated boundary
%-- k = boundary(p,s) , specifies shrink factor s: 0~1.                --
%-- shp = alphaShape(P,a) creates an alpha shape with alpha radius a  --


%--    Loop to find the most compact boundary    -- 

% for s = 0.8:-0.05:0
%     tri  = boundary(p2,s);
%     mesh = surfaceMesh(p2,tri);
%     TF = isWatertight(mesh);    % Check if surface mesh is watertight
%     if  TF
%         figure;
%         trimesh(tri,p2(:,1),p2(:,2),p2(:,3));
%         fprintf('%s\n','Using alpha Shape algorithm generates boundaries');
%         fprintf('%s\n',['Shrink factor=',num2str(s)]);
%         break
%     end
% end


Alpha0        = 1 ;       % initial alpha radius
shp       = alphaShape(p2);
alphaspec = alphaSpectrum(shp) ;%Find the spectrum of critical alpha values that produce unique alpha shapes for the point cloud.

maxa      = max(alphaspec);
da        = min(alphaspec);
for Alpha = Alpha0:da: maxa
    [tri, p2] = boundaryFacets(shp);
    mesh      = surfaceMesh(p2,tri);
    TF        = isWatertight(mesh);    % Check if surface mesh is watertight
    if  TF
        figure;
        plot(shp);
        fprintf('%s\n','Using alpha Shape algorithm generates boundaries');
        fprintf('%s\n',['Alpha radius =',num2str(shp.Alpha)]);
        break
    else
        shp.Alpha = Alpha;
    end
end
 

if  ~TF
    fprintf('%s\n','Using convhull function generates boundaries');
    [tri,av1] = convhull(p2);
    mesh      = surfaceMesh(p2,tri);
    trimesh(tri,p2(:,1),p2(:,2),p2(:,3));
    TF        = isWatertight(mesh);
    if ~TF
        error("surface mesh isn't watertight")
    end
end

% trimesh(tri,xyz(:,1),xyz(:,2),xyz(:,3))
% trisurf(tri,p2(:,1),p2(:,2),p2(:,3),...
%     'FaceColor','cyan','FaceAlpha',1)
% axis equal
% % scatter3(p(:,1),p(:,2),p(:,3))


removeDefects(mesh,"duplicate-vertices");
removeDefects(mesh,"duplicate-faces");
removeDefects(mesh,"unreferenced-vertices");
removeDefects(mesh,"degenerate-faces");
removeDefects(mesh,"nonmanifold-edges");


% Check If Surface Mesh Is Self-Intersecting
if isSelfIntersecting(mesh) 
warning 'At least one of its face intersects another face.' ;
end


subdivide(mesh,"loop",3);%"midpoint-split"
surfaceMeshShow(mesh,Title="Subdivided Mesh",WireFrame=true)
figure
trimesh(mesh.Faces,mesh.Vertices(:,1),mesh.Vertices(:,2),mesh.Vertices(:,3))
colormap("jet")
hold on;scatter3(p2(:,1),p2(:,2),p2(:,3),'filled')

% writeSurfaceMesh(mesh,"ManifoldMesh2.ply")


% Find vertex center of surface mesh
% vertexCenter(mesh);


numIterations = 2;
scaleFactor = [0.2];%[-0.62 0.6]
Mesh = smoothSurfaceMesh(mesh,numIterations,Method="Laplacian",ScaleFactor=scaleFactor);%"Laplacian" "Taubin"
surfaceMeshShow(Mesh,Title="Laplacian")

writeSurfaceMesh(Mesh,"smoothSurfaceMesh.ply");

