clear;clc;close all;
path('toolbox',path);
options.USING_POINT_RING = GS.USING_POINT_RING;
extension='.off';

%% Step 0: read file (point cloud & local feature size if possible), and
% normalize the modle.
% filename = '../data/simplejoint_v4770';% which file we should run on
filename = '../data/ls';
tic
P.filename = [filename extension];% point set
[P.pts,P.faces] = read_mesh(P.filename);
P.npts = size(P.pts,1);
if exist([filename '_fe.txt'],'file') % result of Tamal K Dey's NormFet
    P.radis = load([filename '_fe.txt']);
else
    P.radis = ones(P.npts,1);
end

% P.pts = GS.normalize(P.pts);

            bbox = [min(P.pts(:,1)), min(P.pts(:,2)), min(P.pts(:,3)), max(P.pts(:,1)), max(P.pts(:,2)), max(P.pts(:,3))];
            c = (bbox(4:6)+bbox(1:3))*0.5;            
            P.pts = P.pts - repmat(c, size(P.pts,1), 1);
            s = 1.6 / max(bbox(4:6)-bbox(1:3));% make the bbox's diagnol = 1.6. %1.0, 1.6
            P.pts = P.pts*s;


O_pts=P.pts/s;
O_pts=O_pts+repmat(c, size(P.pts,1), 1);

[P.bbox, P.diameter] = GS.compute_bbox(P.pts);
disp(sprintf('read point set:'));
toc

%% Step 1: build local 1-ring
% build neighborhood, knn?
tic
P.k_knn = GS.compute_k_knn(P.npts);
if options.USING_POINT_RING
    P.rings = compute_point_point_ring(P.pts, P.k_knn, []);
else    
    P.frings = compute_vertex_face_ring(P.faces);
    P.rings = compute_vertex_ring(P.faces, P.frings);
end
disp(sprintf('compute local 1-ring:'));
toc

%% Step 1: Contract point cloud by Laplacian
tic
[P.cpts, t, initWL, WC, sl] = find_center_point(P, options);

O_cpts=P.cpts/s;
O_cpts=O_cpts+repmat(c, size(P.cpts,1), 1);

% center points
Pc_x=O_cpts(:,1);
Pc_y=O_cpts(:,2);
Pc_z=O_cpts(:,3);

% original points
si_x=O_pts(:,1);
si_y=O_pts(:,2);
si_z=O_pts(:,3);


m=1;
Spts_x=zeros(1,300);
Spts_y=zeros(1,300);
Spts_z=zeros(1,300);
for i=1:1926
    D_min=100;
    for j=1:1926
    D=E_dis(Pc_x(i),Pc_y(i),Pc_z(i),si_x(j),si_y(j),si_z(j));
    if (D_min>D)
        D_min=D;
    end
    end
    if (D_min>8)
        Spts_x(m)=Pc_x(i);Spts_y(m)=Pc_y(i);Spts_z(m)=Pc_z(i);
        m=m+1;
    end
end

y=[Spts_x; Spts_y; Spts_z];
fid = fopen('exp.txt', 'w');
fprintf(fid, '%6.2f %6.2f %6.2f\r\n',y);
fclose(fid);


figure(4);
plot3(Pc_x(:),Pc_y(:),Pc_z(:),'k*');

figure(5);
scatter3(O_pts(:,1),O_pts(:,2), O_pts(:,3),10,'b','filled');
hold on;
scatter3(O_cpts(:,1),O_cpts(:,2), O_cpts(:,3),10,'r','filled');

C=load('result1.txt');
figure;
plot3(C(:,1),C(:,2),C(:,3),'k.');

disp(sprintf('Contraction:'));
toc