clear all
close all
clc

%% LIGAR CAMARAS
% LIGAR CAMARA RGB
RGBDisp = videoinput('kinect',1,'RGB_640x480');

% LIGAR CAMARA DEPTH
DEPTHDisp = videoinput('kinect',2,'Depth_640x480');

preview(DEPTHDisp)
pause(20)

RGBFrameFlipped = flip(getsnapshot(RGBDisp),2);
DEPTHFrameFlipped = flip(getsnapshot(DEPTHDisp),2);

%% NUVEM DE PONTOS COMPLETA
% EXTRAIR A NUVEM DE PONTOS
RGBFrame1 = getsnapshot(RGBDisp);
DEPTHFrame1 = getsnapshot(DEPTHDisp);

% figure
% imshow(RGBFrameFlipped)
% figure
% imshow(DEPTHFrameFlipped,[])
% [filledDepth,~] =Kinect_DepthNormalization(DEPTHFrame1);
% figure
% imshow(filledDepth,[])

ptCloudOriginal = pcfromkinect(DEPTHDisp, DEPTHFrame1, RGBFrame1);

%% RODAR A NUVEM DE PONTOS PARA FICAR PARALELA COM EIXO DOS XX
%ANGULO DE ROTAÇAO
theta = -0.04;

% MATRIZ DE ROTAÇAO EM X
Rx = [1 0 0; 0 cos(theta) -sin(theta); 0 sin(theta) cos(theta)];
trans = [0, 0, 0];
tform = rigid3d(Rx,trans);

% NOVA NUVEM DE PONTOS RODADA
ptCloudRot = pctransform(ptCloudOriginal,tform);

ptCloudRot = pcdenoise(ptCloudRot,NumNeighbors=50);

% MOSTAR NUVEM DE PONTOS
PointCloudBrowser = figure("Name",'Point Cloud Browser','NumberTitle','off');
pcshow(ptCloudRot,"VerticalAxis","Y",VerticalAxisDir="Down");
title('Full Point Cloud')
whitebg(gcf,'white')
set(gcf,'Color','w')
xlabel('X(m)')
ylabel('Y(m)')
zlabel('Z(m)')


%% DEFINIR UMA ROI CUBICA PARA ASSUMIR A TARGET AREA
% ESPECIFICAR O STEP SIZE
stepSize = 100;

% AMOSTRA A NUVEM DE PONTOS DE ENTRADA E ARMAZENA AS COORDENADAS DE PONTOS COMO UM OBJETO DE NUVEM DE PONTOS
indices = 1:stepSize:ptCloudRot.Count;
tempPtCloud = select(ptCloudRot,indices);

% REMOVER OS PONTOS INVALIDOS DENTRO DA AMOSTRA ANTERIOR
[tempPtCloud,validIndices] = removeInvalidPoints(tempPtCloud);

% OBTER AS COORDENADAS DE LOCALIZAÇAO 3D DOS PONTOS DA NUVEM DE PONTOS
worldPoints = tempPtCloud.Location;

% CONVERTER E ENCONTRAR AS COORDENADAS DOS PONTOS 2D CORRESPONDENTES AOS PONTOS 3D DA NUVEM
[Y,X] = ind2sub([size(ptCloudRot.Location,1),size(ptCloudRot.Location,2)],indices);
imagePoints = [X(validIndices)' Y(validIndices)'];

% ESTIMATIVA DA MATRIZ DE PROJEÇAO DA CAMARA A PARTIR DA IMAGEM E DAS COORDENADAS DOS PONTOS DE LOCALIZAÇAO 3D
camMatrix = estimateCameraMatrix(imagePoints,worldPoints);

% DEFINIR A ROI CUBICA PARA A TARGET AREA
roi = [-0.46 0.5 -0.6 0.418 1.4 2.28];

% ENCONTRAR OS INDICES DOS PONTOS DA NUVEM QUE ESTAO DENTRO DA ROI CUBICA
indices = findPointsInROI(ptCloudRot,roi,camMatrix);

% ptCloudRot.Location = reshape(ptCloudRot.Location,[480,640,3]);

% SELECIONAR OS PONTOS DA NUVEM QUE SE ENCONTRAM DENTRO DA ROI
ptcloudB = select(ptCloudRot,indices);

% MOSTRAR A NUVEM DE PONTOS COM A TARGER AREA
figure(PointCloudBrowser);
pcshow(ptCloudRot,"VerticalAxis","Y",VerticalAxisDir="Down")
hold on
pcshow(ptcloudB.Location,'r',"VerticalAxis","Y",VerticalAxisDir="Down")
title('Full Point Cloud with Target Area')
xlabel('X(m)')
ylabel('Y(m)')
zlabel('Z(m)')

% DEFINIÇÃO DOS PONTOS DOS VERTICES DO CUBO DA ROI
p1 = [roi(1) roi(4) roi(6)]';
p2 = [roi(1) roi(4) roi(5)]';
p3 = [roi(2) roi(4) roi(5)]';
p4 = [roi(2) roi(4) roi(6)]';
p5 = [roi(1) roi(3) roi(6)]';
p6 = [roi(1) roi(3) roi(5)]';
p7 = [roi(2) roi(3) roi(5)]';
p8 = [roi(2) roi(3) roi(6)]';

% DEFINIR AS FACES DO CUBO ATRAVES DOS VERTICES ANTERIORES
A1=[p1 p2 p3 p4;p4 p3 p7 p8;p7 p8 p5 p6;p7 p6 p2 p3;p6 p5 p1 p2;p8 p5 p1 p4];

% CORES DAS FACES
p=fill3(A1(1,:),A1(2,:),A1(3,:),'--w');
p.FaceAlpha = 0.3;
p=fill3(A1(4,:),A1(5,:),A1(6,:),'--w');
p.FaceAlpha = 0.3;
p=fill3(A1(7,:),A1(8,:),A1(9,:),'--w');
p.FaceAlpha = 0.3;
p=fill3(A1(10,:),A1(11,:),A1(12,:),'--w');
p.FaceAlpha = 0.3;
p=fill3(A1(13,:),A1(14,:),A1(15,:),'--w');
p.FaceAlpha = 0.3;
p=fill3(A1(16,:),A1(17,:),A1(18,:),'--w');
p.FaceAlpha = 0.3;

hold off




%% REMOVER O PLANO DO CHAO E MOSTRAR OBJETO
% DEFINIR A DISTANCIA MAXIMA PONTO-PLANO PARA A MONTAGEM DO PLANO
maxDistance = 0.0005;     % (m)

% DEFINIR O VETOR NORMAL AO PLANO A DEFINIR (NESTE CASO VETOR PREPENDICULAR AO PLANO Y)
referenceVector = [0,1,0];

% DETETAR O PLANO E EXTRAI-LO DA NUVEM DE PONTOS
[~,~,outlierIndices] = pcfitplane(ptcloudB,maxDistance,referenceVector);
remainPtCloud = select(ptcloudB,outlierIndices);

% REMOVER O BARULHO DA NUVEM DE PONTOS (ELIMINAR "IMPUREZAS")
pclimpa = pcdenoise(remainPtCloud);

% MOSTRAR A NUVEM DE PONTOS LIMPA E SEM O PLANO DO CHAO
figure("Name",'Object Cloud Viewer','NumberTitle','off');
pcshow(pclimpa,"VerticalAxis","Y",VerticalAxisDir="Down")
title('Object Cloud Point')
whitebg(gcf,'white')
set(gcf,'Color','w')
xlabel('X(m)')
ylabel('Y(m)')
zlabel('Z(m)')



% %% DESENHAR ROI PARA CORTAR O OBJETO
% % CONVERTER MATRIX FULL EM MATRIX DE 3xCOLUNAS COM AS COORDENADAS X Y E Z RESPECTIVAMENTE
% NovaMatrizNuvemTotal = [reshape(ptCloudRot.Location(:,:,1),[],1) reshape(ptCloudRot.Location(:,:,2),[],1) reshape(ptCloudRot.Location(:,:,3),[],1)];
% 
% % COMPARAR MATRIZ FULL E MATRIX DA TARGET AREA E RETIRAR OS INDICES DOS PONTOS DO SEU INTERIOR
% [RowNova01,~] = find(NovaMatrizNuvemTotal(:,1)<roi(2) & NovaMatrizNuvemTotal(:,1)>roi(1) & NovaMatrizNuvemTotal(:,2)<roi(4) & NovaMatrizNuvemTotal(:,2)>roi(3) & NovaMatrizNuvemTotal(:,3)>roi(5) & NovaMatrizNuvemTotal(:,3)<roi(6));
% 
% % CONVERTER OS INDICES EM LINHA COLUNA DENTRO DA RESOLUÇAO DO FRAME
% [ROWNOVA1,COLNOVA1] = ind2sub([480 640],RowNova01);
% 
% % ENCONTRAR OS MAX E MIN DE COLUNA E LINHA PARA FORMATAR IMAGEM
% maxrow3 = max(ROWNOVA1);
% minrow3 = min(ROWNOVA1);
% maxcol3 = max(COLNOVA1);
% mincol3 = min(COLNOVA1);
% 
% % RETIRAR UM FRAME DA CAMARA RGB E DEPTH
% DEPTHFrameFlipped2 = flip(getsnapshot(DEPTHDisp),2);
% RGBFrameFlipped2 = flip(getsnapshot(RGBDisp),2);
% 
% % "RECORTAR" IMAGEM RGB E DEPTH PELOS MAX/MIN COLUNA/LINHA
% RGB_ROI = RGBFrameFlipped2(minrow3:maxrow3,mincol3:maxcol3,:);
% DEPTH_ROI = DEPTHFrameFlipped2((minrow3:maxrow3)-10,(mincol3:maxcol3)+15,:);
% 
% % ABRIR PLOT PARA CRIAR A ZONA ROI
% figroi = figure("Name",'Region of Interest Selection','NumberTitle','off');
% imshow(RGB_ROI)
% 
% h1 = drawrectangle;
% 
% % ESPERAR PELO USER DEFINIR A ROI PARA SE OBTER A SUA POSIÇAO
% wait(h1)
% 
% % RETIRAR DO OBJETO O INTERIOR DA ROI ANTERIORMENTE DEFINIDA
% % ESPECIFICAR O STEP SIZE
% stepSize3 = 100;
% 
% % AMOSTRA A NUVEM DE PONTOS DE ENTRADA E ARMAZENA AS COORDENADAS DE PONTOS COMO UM OBJETO DE NUVEM DE PONTOS
% indices3 = 1:stepSize3:ptCloudRot.Count;
% tempPtCloud3 = select(ptCloudRot,indices3);
% 
% % REMOVER OS PONTOS INVALIDOS DENTRO DA AMOSTRA ANTERIOR
% [tempPtCloud3,validIndices3] = removeInvalidPoints(tempPtCloud3);
% 
% % OBTER AS COORDENADAS DE LOCALIZAÇAO 3D DOS PONTOS DA NUVEM DE PONTOS
% worldPoints3 = tempPtCloud3.Location;
% 
% % CONVERTER E ENCONTRAR AS COORDENADAS DOS PONTOS 2D CORRESPONDENTES AOS PONTOS 3D DA NUVEM
% [Y,X] = ind2sub([size(ptCloudRot.Location,1),size(ptCloudRot.Location,2)],indices3);
% imagePoints3 = [X(validIndices3)' Y(validIndices3)'];
% 
% % ESTIMATIVA DA MATRIZ DE PROJEÇAO DA CAMARA A PARTIR DA IMAGEM E DAS COORDENADAS DOS PONTOS DE LOCALIZAÇAO 3D
% camMatrix3 = estimateCameraMatrix(imagePoints3,worldPoints3);
% 
% % DEFINIR A ROI CUBICA PARA A TARGET AREA
% distbetween_x = norm(max(pclimpa.Location(:,1))-min(pclimpa.Location(:,1)));
% distbetween_y = norm(max(pclimpa.Location(:,2))-min(pclimpa.Location(:,2)));
% 
% Vertice1xmin = floor(h1.Vertices(1));
% Vertice1xmax = floor(h1.Vertices(4));
% Vertice1ymin = floor(h1.Vertices(5));
% Vertice1ymax = floor(h1.Vertices(7));
% 
% if Vertice1xmin == 0
%     Vertice1xmin = 1;
% end
% 
% if Vertice1xmax == 0
%     Vertice1xmax = 1;
% end
% 
% if Vertice1ymin == 0
%     Vertice1ymin = 1;
% end
% 
% if Vertice1ymax == 0
%     Vertice1ymax = 1;
% end
% 
% P_ROI_Xmin = min(min(pclimpa.Location(:,1)))+((Vertice1xmin)*distbetween_x)/(size(RGB_ROI,2));
% P_ROI_Xmax = min(min(pclimpa.Location(:,1)))+((Vertice1xmax)*distbetween_x)/(size(RGB_ROI,2));
% P_ROI_Ymax = min(min(pclimpa.Location(:,2)))+((Vertice1ymax)*distbetween_y)/(size(RGB_ROI,1));
% P_ROI_Ymin = min(min(pclimpa.Location(:,2)))+((Vertice1ymin)*distbetween_y)/(size(RGB_ROI,1));
% 
% roi3 = [P_ROI_Xmin P_ROI_Xmax P_ROI_Ymin P_ROI_Ymax roi(5) roi(6)];
% 
% % ENCONTRAR OS INDICES DOS PONTOS DA NUVEM QUE ESTAO DENTRO DA ROI CUBICA
% indices3 = findPointsInROI(ptCloudRot,roi3,camMatrix3);
% 
% % SELECIONAR OS PONTOS DA NUVEM QUE SE ENCONTRAM DENTRO DA ROI
% ptcloudB3 = select(ptCloudRot,indices3);
% 
% % FECHAR A JANELA DA ROI DEFINIDA
% close(figure(figroi))
% 
% % MOSTRAR A NUVEM DE PONTOS DA ROI
% figure("Name",'ROI Point Cloud Viewer','NumberTitle','off')
% pcshow(ptcloudB3,"VerticalAxis","Y",VerticalAxisDir="Down")
% title('ROI Point Cloud')
% xlabel('X(m)')
% ylabel('Y(m)')
% zlabel('Z(m)')
% 

%% CONTORNOS

% switch app.DropDownSelectROIForContour.Value
%     case 'Full Object'
% 
%         % DEFINIR O VALOR DO INICIAL DO SLIDER PARA 1
%         raio_strel = 1;
% 
%         % EXTRAIR COORDENADAS X E Y DO OBJETO (RETIRAR Z)
%         xnuvem = pclimpa.Location(:,1);
%         ynuvem = pclimpa.Location(:,2);
% 
%     case 'ROI'
%         % DEFINIR O VALOR DO INICIAL DO SLIDER PARA 1
%         raio_strel = 5;
% 
%         % EXTRAIR COORDENADAS X E Y DO OBJETO (RETIRAR Z)
%         xnuvem = ptcloudB3.Location(:,1);
%         ynuvem = ptcloudB3.Location(:,2);
% end


xnuvem = pclimpa.Location(:,1);
ynuvem = pclimpa.Location(:,2);


% FAZER PLOT DOS PONTOS X E Y DO OBJETO E RODAR O PLOT
figure(100);
plot(xnuvem,ynuvem,'k.')
% set(gcf, 'Position', get(0, 'Screensize'));
set(gca,'ydir','reverse')
axis 'off'


% RETIRAR FRAME DO PLOT ANTERIOR

% FramePlot = frame2im(getframe(figure(100)));
FramePlot = getframe(figure(100));
imgFramePlot = FramePlot.cdata;
figure;
imshow(imgFramePlot)



% BINARIZAR IMG E INVERTER COR
ImgBin = imbinarize(im2gray(imgFramePlot));
InverseImg = uint8(255) - im2uint8(ImgBin);

% PREENCHER ESPAÇOS ABERTOS NO OBJETO
raio_strel = 8;
se = strel('disk',raio_strel);
ImgClosed = imclose(InverseImg,se);
figure
imshow(ImgClosed)



[B,L,N,A] = bwboundaries(ImgClosed); 
% figure
% imshow(ImgClosed)
figure;
% subplot(2,4,3:4)
triplot(DT);
xlabel('x')
ylabel('y')
zlabel('z')
xlim([min(xmeshnova)-5 max(xmeshnova)+5])
ylim([min(ymeshnova)-5 max(ymeshnova)+5])
title('2D Delaunay Mesh')
hold on; 
% Loop through object boundaries  
for k = 1:N
    % Boundary k is the parent of a hole if the k-th column of the adjacency matrix A contains a non-zero element
    if (nnz(A(:,k)) > 0)
        boundary = B{k};
        redplot = plot(boundary(:,2),boundary(:,1),'r','LineWidth',2);
        % Loop through the children of boundary k
        for l = find(A(:,k))'
            boundary = B{l};
           greenplot = plot(boundary(:,2),boundary(:,1),'g','LineWidth',2);
        end
    end
end






figure
axes
plot(redplot.XData,redplot.YData,'r','LineWidth',2)
hold on
plot(greenplot.XData,greenplot.YData,'g','LineWidth',2); 
set(gca,'ydir','reverse')




































object_cont = B{1};
holcont1 = B{2};
holcont2 = B{3};
holcont3 = B{4};

figure
plot(object_cont(:,2),object_cont(:,1))
hold on
plot(holcont1(:,2),holcont1(:,1))
hold on
plot(holcont2(:,2),holcont2(:,1))
hold on
plot(holcont3(:,2),holcont3(:,1))
set(gca,'ydir','reverse')



% [B1,L1,N1] = bwboundaries(ImgClosed);
% figure
% imshow(ImgClosed); hold on;
% for k=1:length(B1)
%    boundary = B1{k};
%    if(k > N1)
%      plot(boundary(:,2), boundary(:,1), 'g','LineWidth',2);
%    else
%      plot(boundary(:,2), boundary(:,1), 'r','LineWidth',2);
%    end
% end



% FAZER CONTORNOS DO OBJETO
ImgEdge = edge(ImgClosed,"canny");
% 
% % ENCONTRAR LINHAS E COLUNAS DOS PONTOS BRANCOS ( AKA CONTORNOS)
% [rowbw1,colbw1] = find(ImgEdge==1);
% 
% % LIMITAR IMAGEM SEGUNDO OS CONTORNOS
% rowmax12 = max(rowbw1)+10;
% rowmin12 = min(rowbw1)-10;
% colmax12 = max(colbw1)+10;
% colmin12 = min(colbw1)-10;
% 
% % NAO DEIXAR QUE A IMAGEM RECORTADA NAO ULTRAPASSE A RESOLUÇAO DA IMAGEM
% if rowmin12 < 0
%     rowmin12 = 1;
% end
% 
% if rowmax12 > size(ImgEdge,1)
%     rowmax12 = size(ImgEdge,1);
% end
% 
% if colmin12 < 0
%     colmin12 = 1;
% end
% 
% if colmax12 > size(ImgEdge,2)
%     colmax12 = size(ImgEdge,2);
% end
% 
% % RECORTAR E MOSTRAR IMAGEM DOS CONTORNOS
% novaimagembw1 = ImgEdge(rowmin12:rowmax12,colmin12:colmax12,:);
figure
imshow(ImgEdge)


%% MESH
xmesh = double(pclimpa.Location(:,1));
ymesh = double(pclimpa.Location(:,2));
zmesh = double(pclimpa.Location(:,3));

% MULTIPLICAR POR 1000 PARA OBTER A MESMA ORDEM DE GRANDEZA
matrixmesh = [xmesh ymesh]*1000;

% CENTRAR A MESH COM A ORIGEM DO REFERENCIAL (0,0)
xmeshnova = matrixmesh(:,1)+abs(min(matrixmesh(:,1)));
ymeshnova = abs(matrixmesh(:,2)-abs(max(matrixmesh(:,2))));

% C = [object_cont(:,2),object_cont(:,1)];

DT = delaunayTriangulation(xmeshnova,ymeshnova);

% FULL 2D DELAUNAY MESH OBJETO
figure;
% subplot(2,4,3:4)
triplot(DT);
xlabel('x')
ylabel('y')
zlabel('z')
xlim([min(xmeshnova)-5 max(xmeshnova)+5])
ylim([min(ymeshnova)-5 max(ymeshnova)+5])
title('2D Delaunay Mesh')
hold on
plot(object_cont(:,2),object_cont(:,1),'r','LineWidth',2)



% INDICES DOS PONTOS DOS CONTORNOS
[y,x] = find(ImgEdge); 

% CENTRAR OS CONTORNOS COM A ORIGEM DO REFERENCIAL (0,0)
x1 = x-abs(min(x));
y1 = y-abs(max(y));

% MAX MALHA
max_x_malha = max(xmeshnova);
max_y_malha = abs(min(ymeshnova));

% MAX CONTORNOS
max_x_contorno = max(x1);
max_y_contorno = abs(min(y1));

% FATOR DE MULTIPLICAÇAO DAS MATRIZES ANTERIORES
% X
fator_x = max_x_malha/max_x_contorno;

% Y
fator_y = max_y_malha/max_y_contorno;

% MATRIZ DOS CONTORNOS DO TAMANHO DA MESH
x2 = x1*fator_x;
y2 = y1*fator_y;

% MOSTRAR CONTORNOS SOBREPOSTOS COM A MESH
hold on
plottedcont = plot(x2,y2,'k.','MarkerSize',5);
hold off



% MALHA CRIADA COM ALPHASHAPE E NAO DELAUNAY
shp = alphaShape(xmeshnova,abs(ymeshnova));
% figure
subplot(2,4,1)
plot(shp,'FaceColor','red')
title('AlphaShape 2D Mesh c/Alpha original')
xlim([min(xmeshnova)-10 max(xmeshnova)+10])
ylim([min(abs(ymeshnova))-10 max(abs(ymeshnova))+10])

% MALHA ALPHASHAPE COM ALPHA ALTERADO
shp.Alpha = 10;
% figure
subplot(2,4,2)
plot(shp,'FaceColor','cyan')
hold on
plot(x2,abs(y2),'k.','MarkerSize',5);
hold off
title('AlphaShape 2D Mesh c/Alpha alterado')
xlim([min(xmeshnova)-10 max(xmeshnova)+10])
ylim([min(abs(ymeshnova))-10 max(abs(ymeshnova))+10])



% MALHA 3D COM ALPHASHAPE
subplot(2,4,5)
shp1 = alphaShape(xmesh,zmesh,ymesh);
plot(shp1,'FaceColor','red')
title('AlphaShape 3D Mesh c/Alpha original')
xlim([min(xmesh) max(xmesh)])
ylim([min(zmesh) max(zmesh)])
zlim([min(ymesh) max(ymesh)])
ax1 = gca;
ax1.ZDir = 'reverse';

% MALHA 3D COM ALPHASHAPE COM ALPHA ALTERADO
subplot(2,4,6)
shp1 = alphaShape(xmesh,zmesh,ymesh);
shp1.Alpha = 0.011;
plot(shp1,'FaceColor','cyan')
title('AlphaShape 3D Mesh c/Alpha alterado')
xlim([min(xmesh) max(xmesh)])
ylim([min(zmesh) max(zmesh)])
zlim([min(ymesh) max(ymesh)])
ax1 = gca;
ax1.ZDir = 'reverse';




% MESH 3D
% figure(500)
subplot(2,4,7:8)
hidden on
trimesh(DT.ConnectivityList,xmesh,zmesh,ymesh)
xlabel('x')
ylabel('z')
zlabel('y')
ax = gca;
ax.ZDir = 'reverse';
title('3D Delaunay Mesh')
xlim([min(xmesh) max(xmesh)])
ylim([min(zmesh) max(zmesh)])
zlim([min(ymesh) max(ymesh)])











generatemesh()


























%% contornos novos


xnuvem1 = pclimpa.Location(:,1);
ynuvem1 = pclimpa.Location(:,2);


% FAZER PLOT DOS PONTOS X E Y DO OBJETO E RODAR O PLOT
figure(1000);
plot(xnuvem1,ynuvem1,'o')
% set(gcf, 'Position', get(0, 'Screensize'));
set(gca,'ydir','reverse')
axis 'off'


% RETIRAR FRAME DO PLOT ANTERIOR

% FramePlot = frame2im(getframe(figure(100)));
FramePlot1 = getframe(figure(1000));
imgFramePlot1 = FramePlot1.cdata;
figure;
imshow(imgFramePlot1)




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% BINARIZAR IMG E INVERTER COR
ImgBin1 = imbinarize(im2gray(imgFramePlot1));
figure
imshow(ImgBin1)

InverseImg1 = uint8(255) - im2uint8(ImgBin1);
figure
imshow(InverseImg1)

% preencher buracos da imagem
imfilled1 = imfill(InverseImg1,'holes');
figure
imshow(imfilled1)


% preencher/dilatar pontos com discos com r raio
se2 = strel("disk",3);
imdilated = imdilate(imfilled1,se2);

figure
imshow(imdilated)



% PREENCHER ESPAÇOS ABERTOS NO OBJETO
raio_strel1 = 15;
se1 = strel('disk',raio_strel1);
ImgClosed1 = imclose(imdilated,se1);
figure
imshow(ImgClosed1)

% FAZER CONTORNOS DO OBJETO

[B1,L1,N1] = bwboundaries(ImgClosed1);

imshow(ImgClosed1); hold on;
for k=1:length(B1)
   boundary = B1{k};
   if(k > N1)
     plot(boundary(:,2), boundary(:,1), 'g','LineWidth',2);
   else
     plot(boundary(:,2), boundary(:,1), 'r','LineWidth',2);
   end
end



figure
contornosplot = plot(boundary(:,2),boundary(:,1));


x_cont = reshape(contornosplot.XData,[],1);
y_cont = reshape(contornosplot.YData,[],1);
contornos = [x_cont,y_cont];

figure
plot(contornos)


ImgEdge1 = edge(ImgClosed1,"canny");%% novos contornos

figure
imshow(ImgEdge1)




%%
























figure
ImgEdgefill = imfill(ImgEdge);




figure
imshow(ImgEdgefill)


[y_imgedge_fill,x_imgedge_fill] = find(ImgEdgefill == 1);




figure
plot(x_imgedge_fill,y_imgedge_fill,'kd')










% tri1 = alphaTriangulation(shp1);







% [xa, ya, index] = nearestConnection(x2, y3);
% figure
% plot(xa,ya,'k')



% [x_c,y_c]=pointsPARAcontour(x2,y2);
% 
% figure
% plot(x_c,y_c)
% set(gca,'ydir','reverse')





% % TENTAR UNIR OS PONTOS DOS CONTORNOS
% [Xout,Yout,varargout]=points2contour(x2,y2,[min(x),max(x)],'ccw');
% 
% figure
% plot(Xout,Yout,'k','MarkerSize',20,'LineWidth',2)
% set(gca,'ydir','reverse')
% 
% hold on
% triplot(DT,xmeshnova,ymeshnova);
% hold off
% 
% C = floor([reshape(Xout,[],1),reshape(Yout,[],1)]);
% 
% x11 = floor(xmeshnova(xmeshnova>0));
% y11 = abs(ymeshnova(1:end-1));
% y12 = y11(y11>0);
% x12 = x11(1:end-14);
% 
% DT1 = delaunayTriangulation(x12,y12,C);

% hold on
% contornoimg = double(novaimagembw1);
 
% [xcont, ycont] = find(novaimagembw1 == 1);
 
% figure
% plot(ycont,xcont,'.')
% set(gca,'ydir','reverse')
 
% novamatrix = [ycont,xcont];
% figure
% plot(novamatrix)
% set(gca,'ydir','reverse')
 
% imshow(novaimagembw1)
% hold off