% ----------------------------------------------------------------------- %
% Author: Mireia Leon Dasi
% European Space Agency - ESAC
% Last modified: 30/07/2020
% Description: matlab script to convert .vtk shape files to .xml files 
% readable by the SBMT software as polygons
% ----------------------------------------------------------------------- %

clc; clear all; close all;

% MODIFY input/output files IDs
Input_file = 'cg-dlr_spg-shap7-v1.0_125Kfacets_regions.vtk';
Output_file_pol = 'cg-dlr_spg-shap7-v1.0_125Kfacets_subregions_pol.txt';
Output_file_pts = 'cg-dlr_spg-shap7-v1.0_125Kfacets_subregions_pts.txt';

% Add savepath
save_path='';

% Read .vtk input file
[pts_cart,face,units] = read_vtk(Input_file);

% Specify number of regions (units) and number of facets  (polygons)
n_units = max(units)+1;
n_pols = size(face,2); 

count = zeros(n_units,1);

% Classify ID of points belonging to each region in region_id
for i = 1:n_units
    for j = 1:n_pols
        if units(j) == i-1
            count(i,1) = count(i,1)+3;
            region_id(i,count(i,1)-2:count(i,1)) = face(:,j)';
        end
    end
end
%%
boundary_id = zeros(n_units,size(region_id,2));

% Classify  ID of points belonging to a boundary in boundary_id
for i = 1:n_units
    M = region_id;
    M(i,:) = [];
    lim = min(find(region_id(i,:)==0))-1;
    count_pts = 0;
    for j = 1:lim
        n_bound = length(find(M ==region_id(i,j)));
        if n_bound >=1
            count_pts = count_pts + 1;
            boundary_id(i,count_pts) = region_id(i,j);
        end
    end
end

% Create a cell to contain the boundary points of each region
region_pts = cell(n_units,1);
region_pts_llr = cell(n_units,1);
Boundary = cell(n_units,1);

%% Convert points from XYZ to lat, long, radius system
pts_llr = cart2llr(pts_cart);

%% Assign the X,Y,Z and Lat,Long,Rad points 
for i = 1:n_units
    lim = min(find(boundary_id(i,:)==0))-1;
    for j = 1:lim
        Pts_cart(:,j) = pts_cart(:,boundary_id(i,j));
       % Pts_llr(:,j) = pts_llr(:,boundary_id(i,j));
    end
    Pts_cart = unique(Pts_cart','rows')';
    %Pts_llr = unique(Pts_llr','rows')';
    region_pts_ordered = order(Pts_cart);
    region_pts_ordered_llr = cart2llr(region_pts_ordered);
    region_pts{i} = region_pts_ordered;
    T = table(region_pts_ordered(1,:)',region_pts_ordered(2,:)',region_pts_ordered(3,:)');
    Output_file = strcat('region_',int2str(i),'.csv');   
    writetable(T,[save_path, Output_file],'Delimiter',' ')
    region_pts_llr{i} = region_pts_ordered_llr;
end

%% Create .txt file for polygons, readable by SBMT
fid = fopen(Output_file_pol,'w');
fprintf(fid,'<?xml version="1.0" encoding="UTF-8" standalone="no"?>\n');
fprintf(fid,'<polygons shapemodel="67P 17442">\n');

formatSpec1 = '    <polygon area="0.0" color="255,0,255" id="1" label="" length="" name="default" vertices="';
formatSpec2 = '"/>\n';
for i = 1:n_units
    [s1,s2] = size(region_pts_llr{i});
    p = reshape(region_pts_llr{i},1,s2*3);
    fprintf(fid,formatSpec1);
    fprintf(fid,'%g ', p);
    fprintf(fid,formatSpec2);
end
fprintf(fid,'</polygons>\n');
fclose(fid)

%% Create .txt file for points, readable by SBMT
id = 0;
%Output_file_pts = 'cg-dlr_spg-shap7-v1.0_50Kfacets_regions_pts.txt';
fid2 = fopen(Output_file_pts,'w');

formatSpec_1 = '%d	default	';
formatSpec_2 = '	NA	NA	NA	NA	0.02	1.0	0.0	255,0,255	""\n';

for i = 1:n_units
    [s1,s2] = size(region_pts_llr{i});
    cart = region_pts{i};
    llr = region_pts_llr{i};
    
    for j = 1:s2
        id = id + 1;

        fprintf(fid2,formatSpec_1,id);
        fprintf(fid2,'%g ', cart(:,j)');
        fprintf(fid2,'%g ', llr(:,j)');
        fprintf(fid2,formatSpec_2);
    end
end
fclose(fid2)


