%% save unique xy points in a file
% This script makes a list of file names of Ottawa upc data, 
% then takes the points from teach file (which are in a 100x100 grid) 
% and puts them into 20x20 cells in 25 different files. 
% It can then read each file, find the unique points, and write out the
% merged set.  The benefit is that the merged set has duplicates removed
% and has better spatial coherence. 
% Points with same xy but differnt z will be consecutive
% This takes 45-90 seconds for each tile, so 25 tiles is probably 30 min.
%% xy origin
% to allow single precision
xorigin =  446300;
yorigin = 5030800;

%% make file names
% tiles = [147:151 166:170 186:190 205:209 225:229]';
% nt = length(tiles);
% f = [repmat('000',nt, 1) num2str(tiles) repmat('.upc',nt,1)];
f = dir('../upc/*.upc');
%% take a funny path through the 5x5 grid when writing out
ipath = [1 2 3 3 3 2 2 1 1 1 1 2 2 3 3 4 5 5 4 4 5 5 4 4 5];
jpath = [1 1 1 2 3 3 2 2 3 4 5 5 4 4 5 5 5 4 4 3 3 2 2 1 1];
%% For each file
for fi =  1:length(f)
    clear xyz
    fname = f(fi).name; % f(fi,:); % f(fi).name;
    disp(fname)
    tic
    fid = fopen(['../upc/' fname]);
    h = readUPChead(fid);
    fout = zeros(5,5);

    while ~feof(fid)  % Read from the file in batches of 100,000 pts
        xyz = readUPCxyz(fid);
        xyz(:,1) = xyz(:,1)*h.xScale+ h.xOffset - xorigin;
        xyz(:,2) = xyz(:,2)*h.yScale+ h.yOffset - yorigin;
        xyz(:,3) = xyz(:,3)*h.zScale+ h.zOffset;
        ji = floor(mod(xyz(:,1:2),100)/20)*[5;1]+6; % 5x5 grid ids

        for i = 1:5
            for j = 1:5
                idx = find(ji==i+j*5); % find those that fall into i,j
                if ~isempty(idx)
                    if fout(i,j) == 0 % open the ouput file if haven't done so.
                        fout(i,j) = fopen(['E:/tmp' num2str(i) '-' num2str(j) '.tmp'],'w');
                        if fout(i,j) < 0
                            error(['Cannot open E:/tmp' num2str(i) '-' num2str(j) '.tmp']);
                        end
                    end
                    fwrite(fout(i,j), xyz(idx,:)', 'double'); % write pts in batch
                end
            end
        end
    end % end while - loop to next batch or done with file
    fclose all;
    toc
    
    % Copy file back and find unique points
    % points with same xy but differnt z will be consecutive
    fid = fopen(['../uniqued/uniqued' fname 'xyz'],'w');
    for k = 1:25
        i = ipath(k);
        j = jpath(k);
        if fout(i,j) > 0
            fin = fopen(['E:/tmp' num2str(i) '-' num2str(j) '.tmp']);
            xyz = fread(fin, [3 Inf], 'double')';
            fclose(fin);
            fwrite(fid,unique(xyz, 'rows')', 'single');
        end
    end
    fclose(fid);
    toc

end