%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2022-01-06 11:30:36
 % @LastEditTime: 2022-01-13 17:39:10

 % cvt energy
%--------------------------
function [E, dE] = cvt_energy(seeds, cpnts, rho, c_edges)

    npnts = size(seeds, 1);
    c_pgon = zeros(npnts, 2); % n*2, [cx, cy]
    rhoP  = zeros(npnts, 1);

    % compute centroid & its rho 
%     figure; polyvec = [];
    parfor i = 1 : npnts
        e = c_edges{i}; % n*4, n*[x1,y1,x2,y2]
        pgon = polyshape(e(:,1), e(:,2));
        [x,y] = centroid(pgon);
        c_pgon(i,:) = [x,y];

        % find the nearest point of current seed
        pi = seeds(i, :);
        dist = vecnorm(cpnts - pi, 2, 2);
        [~, idx] = min(dist);
        rhoP(i) = rho(idx);

%         plot(pgon); hold on;
%         polyvec = [polyvec; pgon];
%         scatter(cpnts(idx, 1), cpnts(idx, 2), "filled"); hold on;
    end
%     plot(polyvec);
%     voronoi(seeds(:,1), seeds(:,2));

    E0  = sum((seeds - c_pgon) .^ 2, 2); % [n,1]
    dE0 = 2 * (seeds - c_pgon);          % [n,2]
    
    E  = sum(E0 .* rhoP);
    dE = dE0 .* rhoP;
    dE(:, 3) = 0;
    dE = reshape(dE', [], 1);
end