%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2022-01-26 20:58:50
 % @LastEditTime: 2022-05-09 20:59:08

 % PARAS:
 %  INPUT: 
 %      indices:      the query indices of seeds
 %  OUTPUT:
 %      c_edges:      voronoi edges for each seed
 %      cneP:         centorid point for each v-cell
 %      vedges:       all edges of the cvt, for the following MMC steps
 %      c_nnbEdgeIDX: the corresponding edge-idx in vedges1 for each seed
 %      c_nnbSeedIDX: neighbor seed-idx for each edge
%--------------------------
function [c_edges, cenP, vedges1, c_nnbEdgeIDX, c_nnbSeedIDX] = intersection_info(...
    seeds, indices, nelx, nely, nelz, microx,optDesign)    

    globalx = nelx * microx;
    globaly = nely * microx;
    globalz = nelz * microx;

    if strcmp(optDesign,'cantilever_hole')
        [c_edges, cenP] = mexVD_file(seeds, indices(:), globalx, globaly, globalz);
    else
        [c_edges, cenP] = mexVD(seeds, indices(:), globalx, globaly, globalz);
    end

    npnts = length(indices);
    vedges0 = [];
    seedidx = [];
    parfor i = 1 : npnts
        e = c_edges{i};
        vedges0 = [vedges0; e];
        seedidx = [seedidx; repmat(indices(i), size(e, 1), 1)];
    end      
    vedges1 = unique_edges(vedges0);

%     figure;
%     x = vedges1(:,[1,4]);
%     y = vedges1(:,[2,5]);
%     z = vedges1(:,[3,6]);
%     plot3(x',y',z','LineWidth',2); hold on;  
%     scatter3(seeds(indices,1), seeds(indices,2), seeds(indices,3), 'filled');

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

        ei = c_edges{i}; % n*6
        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(:,5)-ei(j,5)) < eps & abs(vedges1(:,6)-ei(j,6)) < eps ) | ...
                            (abs(vedges1(:,4)-ei(j,1)) < eps & abs(vedges1(:,5)-ei(j,2)) < eps & abs(vedges1(:,6)-ei(j,3)) < eps & ...
                            abs(vedges1(:,1)-ei(j,4)) < eps & abs(vedges1(:,2)-ei(j,5)) < eps & abs(vedges1(:,3)-ei(j,6)) < eps ) );            
            assert(~isempty(ej));
            if length(ej) > 1
                % find the nearest one
                dist = vecnorm(vedges1(ej, :) - ei(j, :), 2, 2);
                [~, r] = min(dist);
                ej = ej(r);
            end
            edgeid(j) = ej;
        end
        c_nnbEdgeIDX{i} = edgeid;
    end      

    nedge = size(vedges1, 1);
    c_nnbSeedIDX = cell(nedge, 1); % neighbor seed-idx for each edge
    parfor i = 1 : nedge
        ei = vedges1(i, :);
        [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(:,5)-ei(5)) < eps & abs(vedges0(:,6)-ei(6)) < eps ) | ...
                        (abs(vedges0(:,4)-ei(1)) < eps & abs(vedges0(:,5)-ei(2)) < eps & abs(vedges0(:,6)-ei(3)) < eps & ...
                         abs(vedges0(:,1)-ei(4)) < eps & abs(vedges0(:,2)-ei(5)) < eps & abs(vedges0(:,3)-ei(6)) < eps ) );
        c_nnbSeedIDX{i} = seedidx(si);
    end

%     % test
%     figure(4); clf
%     for i = 1 : nedge
%         % voronoi(seeds(:,1), seeds(:,2)); hold on;
% 
%         si = c_nnbSeedIDX{i};
%         scatter3(seeds(si,1), seeds(si,2), seeds(si,3), 'filled'); hold on;
%         
%         plot3(vedges1(i,[1,4]), vedges1(i,[2,5]), vedges1(i,[3,6])); hold on;
%     end
end