%%
clear, close all, clc

%%
[SRTM,R] = readgeoraster('s01_e036_1arc_v3.tif',...
    'OutputType','double');

%%
[LON,LAT] = meshgrid(36:1/360:37,-1:1/360:0);
SRTM = SRTM(1:10:end,1:10:end)/5000;

%%
vertices = [LON(:),LAT(:),SRTM(:)];
faces = delaunay(LON,LAT);
mesh = surfaceMesh(vertices,faces);

%%
zLimits = [min(SRTM(:)),max(SRTM(:))];
demcmap(zLimits)

%%
numFaces = size(faces,1);
cmap = colormap;
numColors = size(cmap,1);

%%
faceColors = zeros(numFaces,3);
for i = 1:numFaces
    vertexIndices = faces(i,:);
    avgZ = mean(SRTM(vertexIndices));
    colorIndex = round((avgZ-zLimits(1))/ ...
        (zLimits(2)-zLimits(1))*(numColors-1))+1;
    colorIndex = min(max(colorIndex,1),numColors);
    faceColors(i,:) = cmap(colorIndex,:);
end

%%
numVertices = size(vertices,1);
vertexColors = zeros(numVertices,3);
vertexCounts = zeros(numVertices,1);

%%
for i = 1:numFaces
    face = faces(i,:);
    color = faceColors(i,:);
    for j = 1:3
        vIdx = face(j);
        vertexColors(vIdx,:) = vertexColors(vIdx,:)+color;
        vertexCounts(vIdx) = vertexCounts(vIdx)+1;
    end
end

%%
vertexCounts(vertexCounts==0) = 1;
vertexColors = vertexColors./vertexCounts;

%%
mesh.VertexColors = vertexColors;
mesh.FaceColors = faceColors;

%%
surfaceMeshShow(mesh)

%%
writeSurfaceMesh(mesh,'srtm_1.ply')
writeSurfaceMesh(mesh,'srtm_2.fbx')

%%
vertices = vertices(~isnan(vertices(:,3)),:);

%%
colorIndices = round(rescale(SRTM(:),1,size(cmap,1)));
colors = cmap(colorIndices,:);

%%
ptCloud = pointCloud(vertices,'color',colors);

%%
pcshow(ptCloud)

%%
pcwrite(ptCloud,'srtm_3.ply','Encoding','binary')

