%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2022-01-13 00:10:54
 % @LastEditTime: 2022-02-24 11:48:06
%--------------------------
function seeds = quadtree_nonregular(rho, cpnts, maxx, minx, maxy, miny, microx, eps, bpoly)
    % 2. compute faces-center points & rho
    % 3. divide into block with min-size/10
    % while any block>min-size:
    %   divide the into 4
    
    min_size = microx / 4;
    max_size = microx / 2;
    
    %% bounding box and divide into block with max-size
    decomp_info = [];
    for i = minx : max_size : maxx-max_size
        x = [i, i + max_size];

        for j = miny : max_size : maxy-max_size
            y = [j, j + max_size];

            [idx, ~] = find(cpnts(:, 1) <= x(2) & cpnts(:, 1) >= x(1) & ...
                            cpnts(:, 2) <= y(2) & cpnts(:, 2) >= y(1));

            if ~isempty(idx)
%                 rho_info = max(rho(idx)) - min(rho(idx));
                rho_info = mean(rho(idx));
                decomp_info = [decomp_info; [rho_info, i,j, max_size]];
            end
        end
    end

    %% quadtree
    flag = 1;
    while (any(decomp_info(:,1) > eps) || any(decomp_info(:,4) > min_size)) && flag
        flag = 0;
        [r, ~] = find(decomp_info(:,1) > eps);
        
        for ri = 1:length(r)
            di = decomp_info(r(ri), :);

            if di(4) > min_size
                flag = 1;
                step = di(4) / 2;
    
                % decomp_info(r(ri), :) = [];
                decomp_info(r(ri), 1) = -1;
    
                for i = di(2) : step : di(2) + di(4) - step
                    x = [i, i+step];
    
                    for j = di(3) : step : di(3) + di(4) - step
                        y = [j, j+step];
    
                        [idx, ~] = find(cpnts(:, 1) <= x(2) & cpnts(:, 1) >= x(1) & ...
                                        cpnts(:, 2) <= y(2) & cpnts(:, 2) >= y(1));
    
                        if ~isempty(idx)
%                             rho_info = max(rho(idx)) - min(rho(idx));
                            rho_info = mean(rho(idx));
                            decomp_info = [decomp_info; [rho_info, i,j, step]];
                        end
                    end
                end
            end
        end
    end

    decomp_info = decomp_info(decomp_info(:,1)>0,:);

    seeds = [decomp_info(:,2) + decomp_info(:,end)/2, decomp_info(:,3) + decomp_info(:,end)/2];    
    figure; voronoi(seeds(:,1), seeds(:,2)); hold on;
%     figure; plot(polyshape(bpoly),'FaceColor','white','EdgeColor','black','LineWidth', 1.5);
%     hold on; scatter(seeds(:,1), seeds(:,2), 8,'filled','k'); axis equal; axis off;

%     % test plot 
%     polyvec = [];
%     for i = 1:size(decomp_info)
%         di = decomp_info(i,:);
%         pgon = [di(2), di(3);
%                 di(2)+di(4), di(3);
%                 di(2)+di(4), di(3)+di(4);
%                 di(2), di(3)+di(4)];
%         polyvec = [polyvec; polyshape(pgon)];
%     end
%     figure; plot(polyvec);    
end