%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2022-01-07 15:42:40
 % @LastEditTime: 2022-01-09 20:31:54

 % PARAS:
 %  INPUT: 
 %      indices:      the query indices of seeds
 %  OUTPUT:
 %      c_edges:      voronoi edges for each seed
 %      vedges:       all edges of the cvt, for the following steps
 %      c_nnbEdgeIDX: the corresponding edge-idx in vedges1 for each seed
 %      nnbSeedIDX:   two neighbor seed-idx for each edge
%--------------------------
function [c_edges, vedges1, c_nnbEdgeIDX, nnbSeedIDX] = intersection_info(seeds, ...
    indices, bpoly)
    
% tic
    c_edges = mexCVT_nonregular(seeds, indices, bpoly);
%     toc
    
    npnts = length(indices);
    vedges0 = [];
    seedidx = [];
    parfor i = 1 : npnts
        e = c_edges{i}; % n*4, n*[x1,y1,x2,y2]
        vedges0 = [vedges0; e];
        seedidx = [seedidx; repmat(indices(i), size(e, 1), 1)];
    end      
    vedges1 = unique_edges(vedges0);

    eps = 1e-5;
    c_nnbEdgeIDX = cell(npnts, 1);
    parfor i = 1 : npnts

        ei = c_edges{i}; % n*4, n*[x1,y1,x2,y2]
        edgeid = zeros(size(ei, 1), 1);

        for j = 1 : size(ei, 1)
            [ej, ~] = find((abs(vedges1(:,1)-ei(j,1)) < eps & abs(vedges1(:,2)-ei(j,2)) < eps & ... 
                            abs(vedges1(:,3)-ei(j,3)) < eps & abs(vedges1(:,4)-ei(j,4)) < eps) | ...
                            (abs(vedges1(:,3)-ei(j,1)) < eps & abs(vedges1(:,4)-ei(j,2)) < eps & ...
                            abs(vedges1(:,1)-ei(j,3)) < eps & abs(vedges1(:,2)-ei(j,4)) < eps) );
            if isempty(ej)
                error("didn't finding edge-ID!");
            end
            edgeid(j) = ej;
        end
        c_nnbEdgeIDX{i} = edgeid;
    end      

    nedge = size(vedges1, 1);
    nnbSeedIDX = zeros(nedge, 2); % two neighbor seed-idx for each edge
    parfor i = 1 : nedge
        ei = vedges1(i, :); % [x1,y1,x2,y2]
        [si, ~] = find((abs(vedges0(:,1)-ei(1)) < eps & abs(vedges0(:,2)-ei(2)) < eps & ... 
                        abs(vedges0(:,3)-ei(3)) < eps & abs(vedges0(:,4)-ei(4)) < eps) | ...
                        (abs(vedges0(:,3)-ei(1)) < eps & abs(vedges0(:,4)-ei(2)) < eps & ...
                        abs(vedges0(:,1)-ei(3)) < eps & abs(vedges0(:,2)-ei(4)) < eps) );

        if length(si) > 1                    
            nnbSeedIDX(i, :) = seedidx(si);
        else
            nnbSeedIDX(i, :) = repmat(seedidx(si), 2, 1);
        end
    end

%     % test
%     for i = 1 : nedge
%         figure(4); clf
%         voronoi(seeds(:,1), seeds(:,2)); hold on;
% 
%         si = nnbSeedIDX(i,:);
%         scatter(seeds(si,1), seeds(si,2), 'filled'); hold on;
%         
%         plot(vedges1(i,[1,3]), vedges1(i,[2,4])); hold on;
%     end
end